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, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
44 #include "gromacs/gmxlib/network.h"
45 #include "gromacs/linearalgebra/nrjac.h"
46 #include "gromacs/math/do_fit.h"
47 #include "gromacs/math/functions.h"
48 #include "gromacs/math/vec.h"
49 #include "gromacs/mdlib/main.h"
50 #include "gromacs/mdtypes/commrec.h"
51 #include "gromacs/mdtypes/fcdata.h"
52 #include "gromacs/mdtypes/inputrec.h"
53 #include "gromacs/mdtypes/mdatom.h"
54 #include "gromacs/mdtypes/state.h"
55 #include "gromacs/pbcutil/ishift.h"
56 #include "gromacs/pbcutil/mshift.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/fatalerror.h"
62 #include "gromacs/utility/pleasecite.h"
63 #include "gromacs/utility/smalloc.h"
65 // TODO This implementation of ensemble orientation restraints is nasty because
66 // a user can't just do multi-sim with single-sim orientation restraints.
68 void init_orires(FILE *fplog
, const gmx_mtop_t
*mtop
,
71 const t_commrec
*cr
, t_oriresdata
*od
,
74 od
->nr
= gmx_mtop_ftype_count(mtop
, F_ORIRES
);
77 /* Not doing orientation restraints */
81 const int numFitParams
= 5;
82 if (od
->nr
<= numFitParams
)
84 gmx_fatal(FARGS
, "The system has %d orientation restraints, but at least %d are required, since there are %d fitting parameters.",
85 od
->nr
, numFitParams
+ 1, numFitParams
);
88 if (ir
->bPeriodicMols
)
90 /* Since we apply fitting, we need to make molecules whole and this
91 * can not be done when periodic molecules are present.
93 gmx_fatal(FARGS
, "Orientation restraints can not be applied when periodic molecules are present in the system");
98 gmx_fatal(FARGS
, "Orientation restraints do not work with MPI parallelization. Choose 1 MPI rank, if possible.");
100 /* Orientation restraints */
107 od
->fc
= ir
->orires_fc
;
114 int *nr_ex
= nullptr;
115 int typeMin
= INT_MAX
;
117 gmx_mtop_ilistloop_t iloop
= gmx_mtop_ilistloop_init(mtop
);
120 while (gmx_mtop_ilistloop_next(iloop
, &il
, &nmol
))
124 gmx_fatal(FARGS
, "Found %d copies of a molecule with orientation restrains while the current code only supports a single copy. If you want to ensemble average, run multiple copies of the system using the multi-sim feature of mdrun.", nmol
);
127 for (int i
= 0; i
< il
[F_ORIRES
].nr
; i
+= 3)
129 int type
= il
[F_ORIRES
].iatoms
[i
];
130 int ex
= mtop
->ffparams
.iparams
[type
].orires
.ex
;
134 for (int j
= od
->nex
; j
< ex
+1; j
++)
140 GMX_ASSERT(nr_ex
, "Check for allocated nr_ex to keep the static analyzer happy");
143 typeMin
= std::min(typeMin
, type
);
144 typeMax
= std::max(typeMax
, type
);
147 /* With domain decomposition we use the type index for indexing in global arrays */
148 GMX_RELEASE_ASSERT(typeMax
- typeMin
+ 1 == od
->nr
, "All orientation restraint parameter entries in the topology should be consecutive");
149 /* Store typeMin so we can index array with the type offset */
150 od
->typeMin
= typeMin
;
152 snew(od
->S
, od
->nex
);
153 /* When not doing time averaging, the instaneous and time averaged data
154 * are indentical and the pointers can point to the same memory.
156 snew(od
->Dinsl
, od
->nr
);
158 const gmx_multisim_t
*ms
= cr
->ms
;
161 snew(od
->Dins
, od
->nr
);
165 od
->Dins
= od
->Dinsl
;
168 if (ir
->orires_tau
== 0)
176 snew(od
->Dtav
, od
->nr
);
177 od
->edt
= std::exp(-ir
->delta_t
/ir
->orires_tau
);
178 od
->edt_1
= 1.0 - od
->edt
;
180 /* Extend the state with the orires history */
181 state
->flags
|= (1<<estORIRE_INITF
);
182 state
->hist
.orire_initf
= 1;
183 state
->flags
|= (1<<estORIRE_DTAV
);
184 state
->hist
.norire_Dtav
= od
->nr
*5;
185 snew(state
->hist
.orire_Dtav
, state
->hist
.norire_Dtav
);
188 snew(od
->oinsl
, od
->nr
);
191 snew(od
->oins
, od
->nr
);
195 od
->oins
= od
->oinsl
;
197 if (ir
->orires_tau
== 0)
203 snew(od
->otav
, od
->nr
);
205 snew(od
->tmpEq
, od
->nex
);
208 for (int i
= 0; i
< mtop
->natoms
; i
++)
210 if (ggrpnr(&mtop
->groups
, egcORFIT
, i
) == 0)
215 snew(od
->mref
, od
->nref
);
216 snew(od
->xref
, od
->nref
);
217 snew(od
->xtmp
, od
->nref
);
219 snew(od
->eig
, od
->nex
*12);
221 /* Determine the reference structure on the master node.
222 * Copy it to the other nodes after checking multi compatibility,
223 * so we are sure the subsystems match before copying.
225 rvec com
= { 0, 0, 0 };
228 gmx_mtop_atomloop_all_t aloop
= gmx_mtop_atomloop_all_init(mtop
);
231 while (gmx_mtop_atomloop_all_next(aloop
, &i
, &atom
))
233 if (mtop
->groups
.grpnr
[egcORFIT
] == nullptr ||
234 mtop
->groups
.grpnr
[egcORFIT
][i
] == 0)
236 /* Not correct for free-energy with changing masses */
237 od
->mref
[j
] = atom
->m
;
238 if (ms
== nullptr || MASTERSIM(ms
))
240 copy_rvec(xref
[i
], od
->xref
[j
]);
241 for (int d
= 0; d
< DIM
; d
++)
243 com
[d
] += od
->mref
[j
]*xref
[i
][d
];
250 svmul(1.0/mtot
, com
, com
);
251 if (ms
== nullptr || MASTERSIM(ms
))
253 for (int j
= 0; j
< od
->nref
; j
++)
255 rvec_dec(od
->xref
[j
], com
);
259 fprintf(fplog
, "Found %d orientation experiments\n", od
->nex
);
260 for (int i
= 0; i
< od
->nex
; i
++)
262 fprintf(fplog
, " experiment %d has %d restraints\n", i
+1, nr_ex
[i
]);
267 fprintf(fplog
, " the fit group consists of %d atoms and has total mass %g\n",
272 fprintf(fplog
, " the orientation restraints are ensemble averaged over %d systems\n", ms
->nsim
);
274 check_multi_int(fplog
, ms
, od
->nr
,
275 "the number of orientation restraints",
277 check_multi_int(fplog
, ms
, od
->nref
,
278 "the number of fit atoms for orientation restraining",
280 check_multi_int(fplog
, ms
, ir
->nsteps
, "nsteps", FALSE
);
281 /* Copy the reference coordinates from the master to the other nodes */
282 gmx_sum_sim(DIM
*od
->nref
, od
->xref
[0], ms
);
285 please_cite(fplog
, "Hess2003");
288 void diagonalize_orires_tensors(t_oriresdata
*od
)
290 if (od
->M
== nullptr)
293 for (int i
= 0; i
< DIM
; i
++)
297 snew(od
->eig_diag
, DIM
);
299 for (int i
= 0; i
< DIM
; i
++)
305 for (int ex
= 0; ex
< od
->nex
; ex
++)
307 /* Rotate the S tensor back to the reference frame */
309 mmul(od
->R
, od
->S
[ex
], TMP
);
310 mtmul(TMP
, od
->R
, S
);
311 for (int i
= 0; i
< DIM
; i
++)
313 for (int j
= 0; j
< DIM
; j
++)
315 od
->M
[i
][j
] = S
[i
][j
];
320 jacobi(od
->M
, DIM
, od
->eig_diag
, od
->v
, &nrot
);
323 for (int i
= 0; i
< DIM
; i
++)
327 for (int i
= 0; i
< DIM
; i
++)
329 for (int j
= i
+1; j
< DIM
; j
++)
331 if (gmx::square(od
->eig_diag
[ord
[j
]]) > gmx::square(od
->eig_diag
[ord
[i
]]))
340 for (int i
= 0; i
< DIM
; i
++)
342 od
->eig
[ex
*12 + i
] = od
->eig_diag
[ord
[i
]];
344 for (int i
= 0; i
< DIM
; i
++)
346 for (int j
= 0; j
< DIM
; j
++)
348 od
->eig
[ex
*12 + 3 + 3*i
+ j
] = od
->v
[j
][ord
[i
]];
354 void print_orires_log(FILE *log
, t_oriresdata
*od
)
358 diagonalize_orires_tensors(od
);
360 for (int ex
= 0; ex
< od
->nex
; ex
++)
362 eig
= od
->eig
+ ex
*12;
363 fprintf(log
, " Orientation experiment %d:\n", ex
+1);
364 fprintf(log
, " order parameter: %g\n", eig
[0]);
365 for (int i
= 0; i
< DIM
; i
++)
367 fprintf(log
, " eig: %6.3f %6.3f %6.3f %6.3f\n",
368 (eig
[0] != 0) ? eig
[i
]/eig
[0] : eig
[i
],
377 real
calc_orires_dev(const gmx_multisim_t
*ms
,
378 int nfa
, const t_iatom forceatoms
[], const t_iparams ip
[],
379 const t_mdatoms
*md
, const rvec x
[], const t_pbc
*pbc
,
380 t_fcdata
*fcd
, history_t
*hist
)
383 real edt
, edt_1
, invn
, pfac
, r2
, invr
, corrfac
, wsv2
, sw
, dev
;
387 rvec
*xref
, *xtmp
, com
, r_unrot
, r
;
390 const real two_thr
= 2.0/3.0;
396 /* This means that this is not the master node */
397 gmx_fatal(FARGS
, "Orientation restraints are only supported on the master rank, use fewer ranks");
400 bTAV
= (od
->edt
!= 0);
411 od
->exp_min_t_tau
= hist
->orire_initf
*edt
;
413 /* Correction factor to correct for the lack of history
416 corrfac
= 1.0/(1.0 - od
->exp_min_t_tau
);
435 for (int i
= 0; i
< md
->nr
; i
++)
437 if (md
->cORF
[i
] == 0)
439 copy_rvec(x
[i
], xtmp
[j
]);
440 mref
[j
] = md
->massT
[i
];
441 for (int d
= 0; d
< DIM
; d
++)
443 com
[d
] += mref
[j
]*xtmp
[j
][d
];
449 svmul(1.0/mtot
, com
, com
);
450 for (int j
= 0; j
< nref
; j
++)
452 rvec_dec(xtmp
[j
], com
);
454 /* Calculate the rotation matrix to rotate x to the reference orientation */
455 calc_fit_R(DIM
, nref
, mref
, xref
, xtmp
, od
->R
);
457 for (int fa
= 0; fa
< nfa
; fa
+= 3)
459 const int type
= forceatoms
[fa
];
460 const int restraintIndex
= type
- od
->typeMin
;
463 pbc_dx_aiuc(pbc
, x
[forceatoms
[fa
+1]], x
[forceatoms
[fa
+2]], r_unrot
);
467 rvec_sub(x
[forceatoms
[fa
+1]], x
[forceatoms
[fa
+2]], r_unrot
);
469 mvmul(od
->R
, r_unrot
, r
);
471 invr
= gmx::invsqrt(r2
);
472 /* Calculate the prefactor for the D tensor, this includes the factor 3! */
473 pfac
= ip
[type
].orires
.c
*invr
*invr
*3;
474 for (int i
= 0; i
< ip
[type
].orires
.power
; i
++)
478 rvec5
&Dinsl
= od
->Dinsl
[restraintIndex
];
479 Dinsl
[0] = pfac
*(2*r
[0]*r
[0] + r
[1]*r
[1] - r2
);
480 Dinsl
[1] = pfac
*(2*r
[0]*r
[1]);
481 Dinsl
[2] = pfac
*(2*r
[0]*r
[2]);
482 Dinsl
[3] = pfac
*(2*r
[1]*r
[1] + r
[0]*r
[0] - r2
);
483 Dinsl
[4] = pfac
*(2*r
[1]*r
[2]);
487 for (int i
= 0; i
< 5; i
++)
489 od
->Dins
[restraintIndex
][i
] = Dinsl
[i
]*invn
;
496 gmx_sum_sim(5*od
->nr
, od
->Dins
[0], ms
);
499 /* Calculate the order tensor S for each experiment via optimization */
500 for (int ex
= 0; ex
< od
->nex
; ex
++)
502 for (int i
= 0; i
< 5; i
++)
504 matEq
[ex
].rhs
[i
] = 0;
505 for (int j
= 0; j
<= i
; j
++)
507 matEq
[ex
].mat
[i
][j
] = 0;
512 for (int fa
= 0; fa
< nfa
; fa
+= 3)
514 const int type
= forceatoms
[fa
];
515 const int restraintIndex
= type
- od
->typeMin
;
516 rvec5
&Dtav
= od
->Dtav
[restraintIndex
];
519 /* Here we update Dtav in t_fcdata using the data in history_t.
520 * Thus the results stay correct when this routine
521 * is called multiple times.
523 for (int i
= 0; i
< 5; i
++)
526 edt
*hist
->orire_Dtav
[restraintIndex
*5 + i
] +
527 edt_1
*od
->Dins
[restraintIndex
][i
];
531 int ex
= ip
[type
].orires
.ex
;
532 real weight
= ip
[type
].orires
.kfac
;
533 /* Calculate the vector rhs and half the matrix T for the 5 equations */
534 for (int i
= 0; i
< 5; i
++)
536 matEq
[ex
].rhs
[i
] += Dtav
[i
]*ip
[type
].orires
.obs
*weight
;
537 for (int j
= 0; j
<= i
; j
++)
539 matEq
[ex
].mat
[i
][j
] += Dtav
[i
]*Dtav
[j
]*weight
;
543 /* Now we have all the data we can calculate S */
544 for (int ex
= 0; ex
< od
->nex
; ex
++)
546 OriresMatEq
&eq
= matEq
[ex
];
547 /* Correct corrfac and copy one half of T to the other half */
548 for (int i
= 0; i
< 5; i
++)
550 eq
.rhs
[i
] *= corrfac
;
551 eq
.mat
[i
][i
] *= gmx::square(corrfac
);
552 for (int j
= 0; j
< i
; j
++)
554 eq
.mat
[i
][j
] *= gmx::square(corrfac
);
555 eq
.mat
[j
][i
] = eq
.mat
[i
][j
];
558 m_inv_gen(&eq
.mat
[0][0], 5, &eq
.mat
[0][0]);
559 /* Calculate the orientation tensor S for this experiment */
560 matrix
&S
= od
->S
[ex
];
566 for (int i
= 0; i
< 5; i
++)
568 S
[0][0] += 1.5*eq
.mat
[0][i
]*eq
.rhs
[i
];
569 S
[0][1] += 1.5*eq
.mat
[1][i
]*eq
.rhs
[i
];
570 S
[0][2] += 1.5*eq
.mat
[2][i
]*eq
.rhs
[i
];
571 S
[1][1] += 1.5*eq
.mat
[3][i
]*eq
.rhs
[i
];
572 S
[1][2] += 1.5*eq
.mat
[4][i
]*eq
.rhs
[i
];
577 S
[2][2] = -S
[0][0] - S
[1][1];
580 const matrix
*S
= od
->S
;
585 for (int fa
= 0; fa
< nfa
; fa
+= 3)
587 const int type
= forceatoms
[fa
];
588 const int restraintIndex
= type
- od
->typeMin
;
589 const int ex
= ip
[type
].orires
.ex
;
591 const rvec5
&Dtav
= od
->Dtav
[restraintIndex
];
592 od
->otav
[restraintIndex
] = two_thr
*
593 corrfac
*(S
[ex
][0][0]*Dtav
[0] + S
[ex
][0][1]*Dtav
[1] +
594 S
[ex
][0][2]*Dtav
[2] + S
[ex
][1][1]*Dtav
[3] +
595 S
[ex
][1][2]*Dtav
[4]);
598 const rvec5
&Dins
= od
->Dins
[restraintIndex
];
599 od
->oins
[restraintIndex
] = two_thr
*
600 (S
[ex
][0][0]*Dins
[0] + S
[ex
][0][1]*Dins
[1] +
601 S
[ex
][0][2]*Dins
[2] + S
[ex
][1][1]*Dins
[3] +
602 S
[ex
][1][2]*Dins
[4]);
606 /* When ensemble averaging is used recalculate the local orientation
607 * for output to the energy file.
609 const rvec5
&Dinsl
= od
->Dinsl
[restraintIndex
];
610 od
->oinsl
[restraintIndex
] = two_thr
*
611 (S
[ex
][0][0]*Dinsl
[0] + S
[ex
][0][1]*Dinsl
[1] +
612 S
[ex
][0][2]*Dinsl
[2] + S
[ex
][1][1]*Dinsl
[3] +
613 S
[ex
][1][2]*Dinsl
[4]);
616 dev
= od
->otav
[restraintIndex
] - ip
[type
].orires
.obs
;
618 wsv2
+= ip
[type
].orires
.kfac
*gmx::square(dev
);
619 sw
+= ip
[type
].orires
.kfac
;
621 od
->rmsdev
= std::sqrt(wsv2
/sw
);
623 /* Rotate the S matrices back, so we get the correct grad(tr(S D)) */
624 for (int ex
= 0; ex
< od
->nex
; ex
++)
627 tmmul(od
->R
, od
->S
[ex
], RS
);
628 mmul(RS
, od
->R
, od
->S
[ex
]);
633 /* Approx. 120*nfa/3 flops */
636 real
orires(int nfa
, const t_iatom forceatoms
[], const t_iparams ip
[],
637 const rvec x
[], rvec4 f
[], rvec fshift
[],
638 const t_pbc
*pbc
, const t_graph
*g
,
639 real gmx_unused lambda
, real gmx_unused
*dvdlambda
,
640 const t_mdatoms gmx_unused
*md
, t_fcdata
*fcd
,
641 int gmx_unused
*global_atom_index
)
643 int ex
, power
, ki
= CENTRAL
;
645 real r2
, invr
, invr2
, fc
, smooth_fc
, dev
, devins
, pfac
;
648 const t_oriresdata
*od
;
656 bTAV
= (od
->edt
!= 0);
661 /* Smoothly switch on the restraining when time averaging is used */
662 smooth_fc
*= (1.0 - od
->exp_min_t_tau
);
665 for (int fa
= 0; fa
< nfa
; fa
+= 3)
667 const int type
= forceatoms
[fa
];
668 const int ai
= forceatoms
[fa
+ 1];
669 const int aj
= forceatoms
[fa
+ 2];
670 const int restraintIndex
= type
- od
->typeMin
;
673 ki
= pbc_dx_aiuc(pbc
, x
[ai
], x
[aj
], r
);
677 rvec_sub(x
[ai
], x
[aj
], r
);
680 invr
= gmx::invsqrt(r2
);
682 ex
= ip
[type
].orires
.ex
;
683 power
= ip
[type
].orires
.power
;
684 fc
= smooth_fc
*ip
[type
].orires
.kfac
;
685 dev
= od
->otav
[restraintIndex
] - ip
[type
].orires
.obs
;
688 * there is no real potential when time averaging is applied
690 vtot
+= 0.5*fc
*gmx::square(dev
);
694 /* Calculate the force as the sqrt of tav times instantaneous */
695 devins
= od
->oins
[restraintIndex
] - ip
[type
].orires
.obs
;
702 dev
= std::sqrt(dev
*devins
);
710 pfac
= fc
*ip
[type
].orires
.c
*invr2
;
711 for (int i
= 0; i
< power
; i
++)
715 mvmul(od
->S
[ex
], r
, Sr
);
716 for (int i
= 0; i
< DIM
; i
++)
719 -pfac
*dev
*(4*Sr
[i
] - 2*(2+power
)*invr2
*iprod(Sr
, r
)*r
[i
]);
724 ivec_sub(SHIFT_IVEC(g
, ai
), SHIFT_IVEC(g
, aj
), dt
);
728 for (int i
= 0; i
< DIM
; i
++)
732 fshift
[ki
][i
] += fij
[i
];
733 fshift
[CENTRAL
][i
] -= fij
[i
];
740 /* Approx. 80*nfa/3 flops */
743 void update_orires_history(t_fcdata
*fcd
, history_t
*hist
)
745 t_oriresdata
*od
= &(fcd
->orires
);
749 /* Copy the new time averages that have been calculated
750 * in calc_orires_dev.
752 hist
->orire_initf
= od
->exp_min_t_tau
;
753 for (int pair
= 0; pair
< od
->nr
; pair
++)
755 for (int i
= 0; i
< 5; i
++)
757 hist
->orire_Dtav
[pair
*5+i
] = od
->Dtav
[pair
][i
];