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
);
90 gmx_fatal(FARGS
, "Orientation restraints do not work with MPI parallelization. Choose 1 MPI rank, if possible.");
92 /* Orientation restraints */
99 od
->fc
= ir
->orires_fc
;
106 int *nr_ex
= nullptr;
107 gmx_mtop_ilistloop_t iloop
= gmx_mtop_ilistloop_init(mtop
);
110 while (gmx_mtop_ilistloop_next(iloop
, &il
, &nmol
))
112 for (int i
= 0; i
< il
[F_ORIRES
].nr
; i
+= 3)
114 int type
= il
[F_ORIRES
].iatoms
[i
];
115 int ex
= mtop
->ffparams
.iparams
[type
].orires
.ex
;
119 for (int j
= od
->nex
; j
< ex
+1; j
++)
128 snew(od
->S
, od
->nex
);
129 /* When not doing time averaging, the instaneous and time averaged data
130 * are indentical and the pointers can point to the same memory.
132 snew(od
->Dinsl
, od
->nr
);
134 const gmx_multisim_t
*ms
= cr
->ms
;
137 snew(od
->Dins
, od
->nr
);
141 od
->Dins
= od
->Dinsl
;
144 if (ir
->orires_tau
== 0)
152 snew(od
->Dtav
, od
->nr
);
153 od
->edt
= std::exp(-ir
->delta_t
/ir
->orires_tau
);
154 od
->edt_1
= 1.0 - od
->edt
;
156 /* Extend the state with the orires history */
157 state
->flags
|= (1<<estORIRE_INITF
);
158 state
->hist
.orire_initf
= 1;
159 state
->flags
|= (1<<estORIRE_DTAV
);
160 state
->hist
.norire_Dtav
= od
->nr
*5;
161 snew(state
->hist
.orire_Dtav
, state
->hist
.norire_Dtav
);
164 snew(od
->oinsl
, od
->nr
);
167 snew(od
->oins
, od
->nr
);
171 od
->oins
= od
->oinsl
;
173 if (ir
->orires_tau
== 0)
179 snew(od
->otav
, od
->nr
);
181 snew(od
->tmp
, od
->nex
);
182 snew(od
->TMP
, od
->nex
);
183 for (int ex
= 0; ex
< od
->nex
; ex
++)
185 snew(od
->TMP
[ex
], 5);
186 for (int i
= 0; i
< 5; i
++)
188 snew(od
->TMP
[ex
][i
], 5);
193 for (int i
= 0; i
< mtop
->natoms
; i
++)
195 if (ggrpnr(&mtop
->groups
, egcORFIT
, i
) == 0)
200 snew(od
->mref
, od
->nref
);
201 snew(od
->xref
, od
->nref
);
202 snew(od
->xtmp
, od
->nref
);
204 snew(od
->eig
, od
->nex
*12);
206 /* Determine the reference structure on the master node.
207 * Copy it to the other nodes after checking multi compatibility,
208 * so we are sure the subsystems match before copying.
210 rvec com
= { 0, 0, 0 };
213 gmx_mtop_atomloop_all_t aloop
= gmx_mtop_atomloop_all_init(mtop
);
216 while (gmx_mtop_atomloop_all_next(aloop
, &i
, &atom
))
218 if (mtop
->groups
.grpnr
[egcORFIT
] == nullptr ||
219 mtop
->groups
.grpnr
[egcORFIT
][i
] == 0)
221 /* Not correct for free-energy with changing masses */
222 od
->mref
[j
] = atom
->m
;
223 if (ms
== nullptr || MASTERSIM(ms
))
225 copy_rvec(xref
[i
], od
->xref
[j
]);
226 for (int d
= 0; d
< DIM
; d
++)
228 com
[d
] += od
->mref
[j
]*xref
[i
][d
];
235 svmul(1.0/mtot
, com
, com
);
236 if (ms
== nullptr || MASTERSIM(ms
))
238 for (int j
= 0; j
< od
->nref
; j
++)
240 rvec_dec(od
->xref
[j
], com
);
244 fprintf(fplog
, "Found %d orientation experiments\n", od
->nex
);
245 for (int i
= 0; i
< od
->nex
; i
++)
247 fprintf(fplog
, " experiment %d has %d restraints\n", i
+1, nr_ex
[i
]);
252 fprintf(fplog
, " the fit group consists of %d atoms and has total mass %g\n",
257 fprintf(fplog
, " the orientation restraints are ensemble averaged over %d systems\n", ms
->nsim
);
259 check_multi_int(fplog
, ms
, od
->nr
,
260 "the number of orientation restraints",
262 check_multi_int(fplog
, ms
, od
->nref
,
263 "the number of fit atoms for orientation restraining",
265 check_multi_int(fplog
, ms
, ir
->nsteps
, "nsteps", FALSE
);
266 /* Copy the reference coordinates from the master to the other nodes */
267 gmx_sum_sim(DIM
*od
->nref
, od
->xref
[0], ms
);
270 please_cite(fplog
, "Hess2003");
273 void diagonalize_orires_tensors(t_oriresdata
*od
)
275 if (od
->M
== nullptr)
278 for (int i
= 0; i
< DIM
; i
++)
282 snew(od
->eig_diag
, DIM
);
284 for (int i
= 0; i
< DIM
; i
++)
290 for (int ex
= 0; ex
< od
->nex
; ex
++)
292 /* Rotate the S tensor back to the reference frame */
294 mmul(od
->R
, od
->S
[ex
], TMP
);
295 mtmul(TMP
, od
->R
, S
);
296 for (int i
= 0; i
< DIM
; i
++)
298 for (int j
= 0; j
< DIM
; j
++)
300 od
->M
[i
][j
] = S
[i
][j
];
305 jacobi(od
->M
, DIM
, od
->eig_diag
, od
->v
, &nrot
);
308 for (int i
= 0; i
< DIM
; i
++)
312 for (int i
= 0; i
< DIM
; i
++)
314 for (int j
= i
+1; j
< DIM
; j
++)
316 if (gmx::square(od
->eig_diag
[ord
[j
]]) > gmx::square(od
->eig_diag
[ord
[i
]]))
325 for (int i
= 0; i
< DIM
; i
++)
327 od
->eig
[ex
*12 + i
] = od
->eig_diag
[ord
[i
]];
329 for (int i
= 0; i
< DIM
; i
++)
331 for (int j
= 0; j
< DIM
; j
++)
333 od
->eig
[ex
*12 + 3 + 3*i
+ j
] = od
->v
[j
][ord
[i
]];
339 void print_orires_log(FILE *log
, t_oriresdata
*od
)
343 diagonalize_orires_tensors(od
);
345 for (int ex
= 0; ex
< od
->nex
; ex
++)
347 eig
= od
->eig
+ ex
*12;
348 fprintf(log
, " Orientation experiment %d:\n", ex
+1);
349 fprintf(log
, " order parameter: %g\n", eig
[0]);
350 for (int i
= 0; i
< DIM
; i
++)
352 fprintf(log
, " eig: %6.3f %6.3f %6.3f %6.3f\n",
353 (eig
[0] != 0) ? eig
[i
]/eig
[0] : eig
[i
],
362 real
calc_orires_dev(const gmx_multisim_t
*ms
,
363 int nfa
, const t_iatom forceatoms
[], const t_iparams ip
[],
364 const t_mdatoms
*md
, const rvec x
[], const t_pbc
*pbc
,
365 t_fcdata
*fcd
, history_t
*hist
)
368 real edt
, edt_1
, invn
, pfac
, r2
, invr
, corrfac
, wsv2
, sw
, dev
;
373 rvec
*xref
, *xtmp
, com
, r_unrot
, r
;
376 const real two_thr
= 2.0/3.0;
382 /* This means that this is not the master node */
383 gmx_fatal(FARGS
, "Orientation restraints are only supported on the master rank, use fewer ranks");
386 bTAV
= (od
->edt
!= 0);
399 od
->exp_min_t_tau
= hist
->orire_initf
*edt
;
401 /* Correction factor to correct for the lack of history
404 corrfac
= 1.0/(1.0 - od
->exp_min_t_tau
);
423 for (int i
= 0; i
< md
->nr
; i
++)
425 if (md
->cORF
[i
] == 0)
427 copy_rvec(x
[i
], xtmp
[j
]);
428 mref
[j
] = md
->massT
[i
];
429 for (int d
= 0; d
< DIM
; d
++)
431 com
[d
] += mref
[j
]*xref
[j
][d
];
437 svmul(1.0/mtot
, com
, com
);
438 for (int j
= 0; j
< nref
; j
++)
440 rvec_dec(xtmp
[j
], com
);
442 /* Calculate the rotation matrix to rotate x to the reference orientation */
443 calc_fit_R(DIM
, nref
, mref
, xref
, xtmp
, R
);
446 /* Index restraint data in order of appearance in forceatoms */
447 int restraintIndex
= 0;
448 for (int fa
= 0; fa
< nfa
; fa
+= 3)
450 int type
= forceatoms
[fa
];
453 pbc_dx_aiuc(pbc
, x
[forceatoms
[fa
+1]], x
[forceatoms
[fa
+2]], r_unrot
);
457 rvec_sub(x
[forceatoms
[fa
+1]], x
[forceatoms
[fa
+2]], r_unrot
);
459 mvmul(R
, r_unrot
, r
);
461 invr
= gmx::invsqrt(r2
);
462 /* Calculate the prefactor for the D tensor, this includes the factor 3! */
463 pfac
= ip
[type
].orires
.c
*invr
*invr
*3;
464 for (int i
= 0; i
< ip
[type
].orires
.power
; i
++)
468 rvec5
&Dinsl
= od
->Dinsl
[restraintIndex
];
469 Dinsl
[0] = pfac
*(2*r
[0]*r
[0] + r
[1]*r
[1] - r2
);
470 Dinsl
[1] = pfac
*(2*r
[0]*r
[1]);
471 Dinsl
[2] = pfac
*(2*r
[0]*r
[2]);
472 Dinsl
[3] = pfac
*(2*r
[1]*r
[1] + r
[0]*r
[0] - r2
);
473 Dinsl
[4] = pfac
*(2*r
[1]*r
[2]);
477 for (int i
= 0; i
< 5; i
++)
479 od
->Dins
[restraintIndex
][i
] = Dinsl
[i
]*invn
;
488 gmx_sum_sim(5*od
->nr
, od
->Dins
[0], ms
);
491 /* Calculate the order tensor S for each experiment via optimization */
492 for (int ex
= 0; ex
< od
->nex
; ex
++)
494 for (int i
= 0; i
< 5; i
++)
497 for (int j
= 0; j
<= i
; j
++)
504 /* Index restraint data in order of appearance in forceatoms */
506 for (int fa
= 0; fa
< nfa
; fa
+= 3)
508 rvec5
&Dtav
= od
->Dtav
[restraintIndex
];
511 /* Here we update Dtav in t_fcdata using the data in history_t.
512 * Thus the results stay correct when this routine
513 * is called multiple times.
515 for (int i
= 0; i
< 5; i
++)
518 edt
*hist
->orire_Dtav
[restraintIndex
*5 + i
] +
519 edt_1
*od
->Dins
[restraintIndex
][i
];
523 int type
= forceatoms
[fa
];
524 int ex
= ip
[type
].orires
.ex
;
525 real weight
= ip
[type
].orires
.kfac
;
526 /* Calculate the vector rhs and half the matrix T for the 5 equations */
527 for (int i
= 0; i
< 5; i
++)
529 rhs
[ex
][i
] += Dtav
[i
]*ip
[type
].orires
.obs
*weight
;
530 for (int j
= 0; j
<= i
; j
++)
532 T
[ex
][i
][j
] += Dtav
[i
]*Dtav
[j
]*weight
;
537 /* Now we have all the data we can calculate S */
538 for (int ex
= 0; ex
< od
->nex
; ex
++)
540 /* Correct corrfac and copy one half of T to the other half */
541 for (int i
= 0; i
< 5; i
++)
543 rhs
[ex
][i
] *= corrfac
;
544 T
[ex
][i
][i
] *= gmx::square(corrfac
);
545 for (int j
= 0; j
< i
; j
++)
547 T
[ex
][i
][j
] *= gmx::square(corrfac
);
548 T
[ex
][j
][i
] = T
[ex
][i
][j
];
551 m_inv_gen(T
[ex
], 5, T
[ex
]);
552 /* Calculate the orientation tensor S for this experiment */
558 for (int i
= 0; i
< 5; i
++)
560 S
[ex
][0][0] += 1.5*T
[ex
][0][i
]*rhs
[ex
][i
];
561 S
[ex
][0][1] += 1.5*T
[ex
][1][i
]*rhs
[ex
][i
];
562 S
[ex
][0][2] += 1.5*T
[ex
][2][i
]*rhs
[ex
][i
];
563 S
[ex
][1][1] += 1.5*T
[ex
][3][i
]*rhs
[ex
][i
];
564 S
[ex
][1][2] += 1.5*T
[ex
][4][i
]*rhs
[ex
][i
];
566 S
[ex
][1][0] = S
[ex
][0][1];
567 S
[ex
][2][0] = S
[ex
][0][2];
568 S
[ex
][2][1] = S
[ex
][1][2];
569 S
[ex
][2][2] = -S
[ex
][0][0] - S
[ex
][1][1];
575 /* Index restraint data in order of appearance in forceatoms */
577 for (int fa
= 0; fa
< nfa
; fa
+= 3)
579 int type
= forceatoms
[fa
];
580 int ex
= ip
[type
].orires
.ex
;
582 const rvec5
&Dtav
= od
->Dtav
[restraintIndex
];
583 od
->otav
[restraintIndex
] = two_thr
*
584 corrfac
*(S
[ex
][0][0]*Dtav
[0] + S
[ex
][0][1]*Dtav
[1] +
585 S
[ex
][0][2]*Dtav
[2] + S
[ex
][1][1]*Dtav
[3] +
586 S
[ex
][1][2]*Dtav
[4]);
589 const rvec5
&Dins
= od
->Dins
[restraintIndex
];
590 od
->oins
[restraintIndex
] = two_thr
*
591 (S
[ex
][0][0]*Dins
[0] + S
[ex
][0][1]*Dins
[1] +
592 S
[ex
][0][2]*Dins
[2] + S
[ex
][1][1]*Dins
[3] +
593 S
[ex
][1][2]*Dins
[4]);
597 /* When ensemble averaging is used recalculate the local orientation
598 * for output to the energy file.
600 const rvec5
&Dinsl
= od
->Dinsl
[restraintIndex
];
601 od
->oinsl
[restraintIndex
] = two_thr
*
602 (S
[ex
][0][0]*Dinsl
[0] + S
[ex
][0][1]*Dinsl
[1] +
603 S
[ex
][0][2]*Dinsl
[2] + S
[ex
][1][1]*Dinsl
[3] +
604 S
[ex
][1][2]*Dinsl
[4]);
607 dev
= od
->otav
[restraintIndex
] - ip
[type
].orires
.obs
;
609 wsv2
+= ip
[type
].orires
.kfac
*gmx::square(dev
);
610 sw
+= ip
[type
].orires
.kfac
;
614 od
->rmsdev
= std::sqrt(wsv2
/sw
);
616 /* Rotate the S matrices back, so we get the correct grad(tr(S D)) */
617 for (int ex
= 0; ex
< od
->nex
; ex
++)
619 tmmul(R
, S
[ex
], TMP
);
625 /* Approx. 120*nfa/3 flops */
628 real
orires(int nfa
, const t_iatom forceatoms
[], const t_iparams ip
[],
629 const rvec x
[], rvec4 f
[], rvec fshift
[],
630 const t_pbc
*pbc
, const t_graph
*g
,
631 real gmx_unused lambda
, real gmx_unused
*dvdlambda
,
632 const t_mdatoms gmx_unused
*md
, t_fcdata
*fcd
,
633 int gmx_unused
*global_atom_index
)
635 int ex
, power
, ki
= CENTRAL
;
637 real r2
, invr
, invr2
, fc
, smooth_fc
, dev
, devins
, pfac
;
640 const t_oriresdata
*od
;
648 bTAV
= (od
->edt
!= 0);
653 /* Smoothly switch on the restraining when time averaging is used */
654 smooth_fc
*= (1.0 - od
->exp_min_t_tau
);
657 /* Index restraint data in order of appearance in forceatoms */
658 int restraintIndex
= 0;
659 for (int fa
= 0; fa
< nfa
; fa
+= 3)
661 int type
= forceatoms
[fa
];
662 int ai
= forceatoms
[fa
+ 1];
663 int aj
= forceatoms
[fa
+ 2];
666 ki
= pbc_dx_aiuc(pbc
, x
[ai
], x
[aj
], r
);
670 rvec_sub(x
[ai
], x
[aj
], r
);
673 invr
= gmx::invsqrt(r2
);
675 ex
= ip
[type
].orires
.ex
;
676 power
= ip
[type
].orires
.power
;
677 fc
= smooth_fc
*ip
[type
].orires
.kfac
;
678 dev
= od
->otav
[restraintIndex
] - ip
[type
].orires
.obs
;
681 * there is no real potential when time averaging is applied
683 vtot
+= 0.5*fc
*gmx::square(dev
);
687 /* Calculate the force as the sqrt of tav times instantaneous */
688 devins
= od
->oins
[restraintIndex
] - ip
[type
].orires
.obs
;
695 dev
= std::sqrt(dev
*devins
);
703 pfac
= fc
*ip
[type
].orires
.c
*invr2
;
704 for (int i
= 0; i
< power
; i
++)
708 mvmul(od
->S
[ex
], r
, Sr
);
709 for (int i
= 0; i
< DIM
; i
++)
712 -pfac
*dev
*(4*Sr
[i
] - 2*(2+power
)*invr2
*iprod(Sr
, r
)*r
[i
]);
717 ivec_sub(SHIFT_IVEC(g
, ai
), SHIFT_IVEC(g
, aj
), dt
);
721 for (int i
= 0; i
< DIM
; i
++)
725 fshift
[ki
][i
] += fij
[i
];
726 fshift
[CENTRAL
][i
] -= fij
[i
];
735 /* Approx. 80*nfa/3 flops */
738 void update_orires_history(t_fcdata
*fcd
, history_t
*hist
)
740 t_oriresdata
*od
= &(fcd
->orires
);
744 /* Copy the new time averages that have been calculated
745 * in calc_orires_dev.
747 hist
->orire_initf
= od
->exp_min_t_tau
;
748 for (int pair
= 0; pair
< od
->nr
; pair
++)
750 for (int i
= 0; i
< 5; i
++)
752 hist
->orire_Dtav
[pair
*5+i
] = od
->Dtav
[pair
][i
];