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.
43 #include "gromacs/gmxlib/network.h"
44 #include "gromacs/linearalgebra/nrjac.h"
45 #include "gromacs/math/do_fit.h"
46 #include "gromacs/math/functions.h"
47 #include "gromacs/math/vec.h"
48 #include "gromacs/mdlib/main.h"
49 #include "gromacs/mdtypes/commrec.h"
50 #include "gromacs/mdtypes/fcdata.h"
51 #include "gromacs/mdtypes/inputrec.h"
52 #include "gromacs/mdtypes/mdatom.h"
53 #include "gromacs/mdtypes/state.h"
54 #include "gromacs/pbcutil/ishift.h"
55 #include "gromacs/pbcutil/mshift.h"
56 #include "gromacs/pbcutil/pbc.h"
57 #include "gromacs/topology/ifunc.h"
58 #include "gromacs/topology/mtop_util.h"
59 #include "gromacs/topology/topology.h"
60 #include "gromacs/utility/fatalerror.h"
61 #include "gromacs/utility/pleasecite.h"
62 #include "gromacs/utility/smalloc.h"
64 // TODO This implementation of ensemble orientation restraints is nasty because
65 // a user can't just do multi-sim with single-sim orientation restraints.
67 void init_orires(FILE *fplog
, const gmx_mtop_t
*mtop
,
70 const t_commrec
*cr
, t_oriresdata
*od
,
73 int i
, j
, d
, ex
, nmol
, *nr_ex
;
76 gmx_mtop_ilistloop_t iloop
;
78 gmx_mtop_atomloop_all_t aloop
;
80 const gmx_multisim_t
*ms
;
82 od
->nr
= gmx_mtop_ftype_count(mtop
, F_ORIRES
);
85 /* Not doing orientation restraints */
89 const int numFitParams
= 5;
90 if (od
->nr
<= numFitParams
)
92 gmx_fatal(FARGS
, "The system has %d orientation restraints, but at least %d are required, since there are %d fitting parameters.",
93 od
->nr
, numFitParams
+ 1, numFitParams
);
98 gmx_fatal(FARGS
, "Orientation restraints do not work with MPI parallelization. Choose 1 MPI rank, if possible.");
100 /* Orientation restraints */
108 od
->fc
= ir
->orires_fc
;
117 iloop
= gmx_mtop_ilistloop_init(mtop
);
118 while (gmx_mtop_ilistloop_next(iloop
, &il
, &nmol
))
120 for (i
= 0; i
< il
[F_ORIRES
].nr
; i
+= 3)
122 ex
= mtop
->ffparams
.iparams
[il
[F_ORIRES
].iatoms
[i
]].orires
.ex
;
126 for (j
= od
->nex
; j
< ex
+1; j
++)
135 snew(od
->S
, od
->nex
);
136 /* When not doing time averaging, the instaneous and time averaged data
137 * are indentical and the pointers can point to the same memory.
139 snew(od
->Dinsl
, od
->nr
);
142 snew(od
->Dins
, od
->nr
);
146 od
->Dins
= od
->Dinsl
;
149 if (ir
->orires_tau
== 0)
157 snew(od
->Dtav
, od
->nr
);
158 od
->edt
= std::exp(-ir
->delta_t
/ir
->orires_tau
);
159 od
->edt_1
= 1.0 - od
->edt
;
161 /* Extend the state with the orires history */
162 state
->flags
|= (1<<estORIRE_INITF
);
163 state
->hist
.orire_initf
= 1;
164 state
->flags
|= (1<<estORIRE_DTAV
);
165 state
->hist
.norire_Dtav
= od
->nr
*5;
166 snew(state
->hist
.orire_Dtav
, state
->hist
.norire_Dtav
);
169 snew(od
->oinsl
, od
->nr
);
172 snew(od
->oins
, od
->nr
);
176 od
->oins
= od
->oinsl
;
178 if (ir
->orires_tau
== 0)
184 snew(od
->otav
, od
->nr
);
186 snew(od
->tmp
, od
->nex
);
187 snew(od
->TMP
, od
->nex
);
188 for (ex
= 0; ex
< od
->nex
; ex
++)
190 snew(od
->TMP
[ex
], 5);
191 for (i
= 0; i
< 5; i
++)
193 snew(od
->TMP
[ex
][i
], 5);
198 for (i
= 0; i
< mtop
->natoms
; i
++)
200 if (ggrpnr(&mtop
->groups
, egcORFIT
, i
) == 0)
205 snew(od
->mref
, od
->nref
);
206 snew(od
->xref
, od
->nref
);
207 snew(od
->xtmp
, od
->nref
);
209 snew(od
->eig
, od
->nex
*12);
211 /* Determine the reference structure on the master node.
212 * Copy it to the other nodes after checking multi compatibility,
213 * so we are sure the subsystems match before copying.
218 aloop
= gmx_mtop_atomloop_all_init(mtop
);
219 while (gmx_mtop_atomloop_all_next(aloop
, &i
, &atom
))
221 if (mtop
->groups
.grpnr
[egcORFIT
] == NULL
||
222 mtop
->groups
.grpnr
[egcORFIT
][i
] == 0)
224 /* Not correct for free-energy with changing masses */
225 od
->mref
[j
] = atom
->m
;
226 if (ms
== NULL
|| MASTERSIM(ms
))
228 copy_rvec(xref
[i
], od
->xref
[j
]);
229 for (d
= 0; d
< DIM
; d
++)
231 com
[d
] += od
->mref
[j
]*xref
[i
][d
];
238 svmul(1.0/mtot
, com
, com
);
239 if (ms
== NULL
|| MASTERSIM(ms
))
241 for (j
= 0; j
< od
->nref
; j
++)
243 rvec_dec(od
->xref
[j
], com
);
247 fprintf(fplog
, "Found %d orientation experiments\n", od
->nex
);
248 for (i
= 0; i
< od
->nex
; i
++)
250 fprintf(fplog
, " experiment %d has %d restraints\n", i
+1, nr_ex
[i
]);
255 fprintf(fplog
, " the fit group consists of %d atoms and has total mass %g\n",
260 fprintf(fplog
, " the orientation restraints are ensemble averaged over %d systems\n", ms
->nsim
);
262 check_multi_int(fplog
, ms
, od
->nr
,
263 "the number of orientation restraints",
265 check_multi_int(fplog
, ms
, od
->nref
,
266 "the number of fit atoms for orientation restraining",
268 check_multi_int(fplog
, ms
, ir
->nsteps
, "nsteps", FALSE
);
269 /* Copy the reference coordinates from the master to the other nodes */
270 gmx_sum_sim(DIM
*od
->nref
, od
->xref
[0], ms
);
273 please_cite(fplog
, "Hess2003");
276 void diagonalize_orires_tensors(t_oriresdata
*od
)
278 int ex
, i
, j
, nrot
, ord
[DIM
], t
;
284 for (i
= 0; i
< DIM
; i
++)
288 snew(od
->eig_diag
, DIM
);
290 for (i
= 0; i
< DIM
; i
++)
296 for (ex
= 0; ex
< od
->nex
; ex
++)
298 /* Rotate the S tensor back to the reference frame */
299 mmul(od
->R
, od
->S
[ex
], TMP
);
300 mtmul(TMP
, od
->R
, S
);
301 for (i
= 0; i
< DIM
; i
++)
303 for (j
= 0; j
< DIM
; j
++)
305 od
->M
[i
][j
] = S
[i
][j
];
309 jacobi(od
->M
, DIM
, od
->eig_diag
, od
->v
, &nrot
);
311 for (i
= 0; i
< DIM
; i
++)
315 for (i
= 0; i
< DIM
; i
++)
317 for (j
= i
+1; j
< DIM
; j
++)
319 if (gmx::square(od
->eig_diag
[ord
[j
]]) > gmx::square(od
->eig_diag
[ord
[i
]]))
328 for (i
= 0; i
< DIM
; i
++)
330 od
->eig
[ex
*12 + i
] = od
->eig_diag
[ord
[i
]];
332 for (i
= 0; i
< DIM
; i
++)
334 for (j
= 0; j
< DIM
; j
++)
336 od
->eig
[ex
*12 + 3 + 3*i
+ j
] = od
->v
[j
][ord
[i
]];
342 void print_orires_log(FILE *log
, t_oriresdata
*od
)
347 diagonalize_orires_tensors(od
);
349 for (ex
= 0; ex
< od
->nex
; ex
++)
351 eig
= od
->eig
+ ex
*12;
352 fprintf(log
, " Orientation experiment %d:\n", ex
+1);
353 fprintf(log
, " order parameter: %g\n", eig
[0]);
354 for (i
= 0; i
< DIM
; i
++)
356 fprintf(log
, " eig: %6.3f %6.3f %6.3f %6.3f\n",
357 (eig
[0] != 0) ? eig
[i
]/eig
[0] : eig
[i
],
366 real
calc_orires_dev(const gmx_multisim_t
*ms
,
367 int nfa
, const t_iatom forceatoms
[], const t_iparams ip
[],
368 const t_mdatoms
*md
, const rvec x
[], const t_pbc
*pbc
,
369 t_fcdata
*fcd
, history_t
*hist
)
371 int fa
, d
, i
, j
, type
, ex
, nref
;
372 real edt
, edt_1
, invn
, pfac
, r2
, invr
, corrfac
, weight
, wsv2
, sw
, dev
;
374 rvec5
*Dinsl
, *Dins
, *Dtav
, *rhs
;
377 rvec
*xref
, *xtmp
, com
, r_unrot
, r
;
380 const real two_thr
= 2.0/3.0;
386 /* This means that this is not the master node */
387 gmx_fatal(FARGS
, "Orientation restraints are only supported on the master rank, use fewer ranks");
390 bTAV
= (od
->edt
!= 0);
406 od
->exp_min_t_tau
= hist
->orire_initf
*edt
;
408 /* Correction factor to correct for the lack of history
411 corrfac
= 1.0/(1.0 - od
->exp_min_t_tau
);
430 for (i
= 0; i
< md
->nr
; i
++)
432 if (md
->cORF
[i
] == 0)
434 copy_rvec(x
[i
], xtmp
[j
]);
435 mref
[j
] = md
->massT
[i
];
436 for (d
= 0; d
< DIM
; d
++)
438 com
[d
] += mref
[j
]*xtmp
[j
][d
];
444 svmul(1.0/mtot
, com
, com
);
445 for (j
= 0; j
< nref
; j
++)
447 rvec_dec(xtmp
[j
], com
);
449 /* Calculate the rotation matrix to rotate x to the reference orientation */
450 calc_fit_R(DIM
, nref
, mref
, xref
, xtmp
, R
);
454 for (fa
= 0; fa
< nfa
; fa
+= 3)
456 type
= forceatoms
[fa
];
459 pbc_dx_aiuc(pbc
, x
[forceatoms
[fa
+1]], x
[forceatoms
[fa
+2]], r_unrot
);
463 rvec_sub(x
[forceatoms
[fa
+1]], x
[forceatoms
[fa
+2]], r_unrot
);
465 mvmul(R
, r_unrot
, r
);
467 invr
= gmx::invsqrt(r2
);
468 /* Calculate the prefactor for the D tensor, this includes the factor 3! */
469 pfac
= ip
[type
].orires
.c
*invr
*invr
*3;
470 for (i
= 0; i
< ip
[type
].orires
.power
; i
++)
474 Dinsl
[d
][0] = pfac
*(2*r
[0]*r
[0] + r
[1]*r
[1] - r2
);
475 Dinsl
[d
][1] = pfac
*(2*r
[0]*r
[1]);
476 Dinsl
[d
][2] = pfac
*(2*r
[0]*r
[2]);
477 Dinsl
[d
][3] = pfac
*(2*r
[1]*r
[1] + r
[0]*r
[0] - r2
);
478 Dinsl
[d
][4] = pfac
*(2*r
[1]*r
[2]);
482 for (i
= 0; i
< 5; i
++)
484 Dins
[d
][i
] = Dinsl
[d
][i
]*invn
;
493 gmx_sum_sim(5*od
->nr
, Dins
[0], ms
);
496 /* Calculate the order tensor S for each experiment via optimization */
497 for (ex
= 0; ex
< od
->nex
; ex
++)
499 for (i
= 0; i
< 5; i
++)
502 for (j
= 0; j
<= i
; j
++)
509 for (fa
= 0; fa
< nfa
; fa
+= 3)
513 /* Here we update Dtav in t_fcdata using the data in history_t.
514 * Thus the results stay correct when this routine
515 * is called multiple times.
517 for (i
= 0; i
< 5; i
++)
519 Dtav
[d
][i
] = edt
*hist
->orire_Dtav
[d
*5+i
] + edt_1
*Dins
[d
][i
];
523 type
= forceatoms
[fa
];
524 ex
= ip
[type
].orires
.ex
;
525 weight
= ip
[type
].orires
.kfac
;
526 /* Calculate the vector rhs and half the matrix T for the 5 equations */
527 for (i
= 0; i
< 5; i
++)
529 rhs
[ex
][i
] += Dtav
[d
][i
]*ip
[type
].orires
.obs
*weight
;
530 for (j
= 0; j
<= i
; j
++)
532 T
[ex
][i
][j
] += Dtav
[d
][i
]*Dtav
[d
][j
]*weight
;
537 /* Now we have all the data we can calculate S */
538 for (ex
= 0; ex
< od
->nex
; ex
++)
540 /* Correct corrfac and copy one half of T to the other half */
541 for (i
= 0; i
< 5; i
++)
543 rhs
[ex
][i
] *= corrfac
;
544 T
[ex
][i
][i
] *= gmx::square(corrfac
);
545 for (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 (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];
576 for (fa
= 0; fa
< nfa
; fa
+= 3)
578 type
= forceatoms
[fa
];
579 ex
= ip
[type
].orires
.ex
;
581 od
->otav
[d
] = two_thr
*
582 corrfac
*(S
[ex
][0][0]*Dtav
[d
][0] + S
[ex
][0][1]*Dtav
[d
][1] +
583 S
[ex
][0][2]*Dtav
[d
][2] + S
[ex
][1][1]*Dtav
[d
][3] +
584 S
[ex
][1][2]*Dtav
[d
][4]);
587 od
->oins
[d
] = two_thr
*(S
[ex
][0][0]*Dins
[d
][0] + S
[ex
][0][1]*Dins
[d
][1] +
588 S
[ex
][0][2]*Dins
[d
][2] + S
[ex
][1][1]*Dins
[d
][3] +
589 S
[ex
][1][2]*Dins
[d
][4]);
593 /* When ensemble averaging is used recalculate the local orientation
594 * for output to the energy file.
596 od
->oinsl
[d
] = two_thr
*
597 (S
[ex
][0][0]*Dinsl
[d
][0] + S
[ex
][0][1]*Dinsl
[d
][1] +
598 S
[ex
][0][2]*Dinsl
[d
][2] + S
[ex
][1][1]*Dinsl
[d
][3] +
599 S
[ex
][1][2]*Dinsl
[d
][4]);
602 dev
= od
->otav
[d
] - ip
[type
].orires
.obs
;
604 wsv2
+= ip
[type
].orires
.kfac
*gmx::square(dev
);
605 sw
+= ip
[type
].orires
.kfac
;
609 od
->rmsdev
= std::sqrt(wsv2
/sw
);
611 /* Rotate the S matrices back, so we get the correct grad(tr(S D)) */
612 for (ex
= 0; ex
< od
->nex
; ex
++)
614 tmmul(R
, S
[ex
], TMP
);
620 /* Approx. 120*nfa/3 flops */
623 real
orires(int nfa
, const t_iatom forceatoms
[], const t_iparams ip
[],
624 const rvec x
[], rvec4 f
[], rvec fshift
[],
625 const t_pbc
*pbc
, const t_graph
*g
,
626 real gmx_unused lambda
, real gmx_unused
*dvdlambda
,
627 const t_mdatoms gmx_unused
*md
, t_fcdata
*fcd
,
628 int gmx_unused
*global_atom_index
)
631 int fa
, d
, i
, type
, ex
, power
, ki
= CENTRAL
;
633 real r2
, invr
, invr2
, fc
, smooth_fc
, dev
, devins
, pfac
;
636 const t_oriresdata
*od
;
644 bTAV
= (od
->edt
!= 0);
649 /* Smoothly switch on the restraining when time averaging is used */
650 smooth_fc
*= (1.0 - od
->exp_min_t_tau
);
654 for (fa
= 0; fa
< nfa
; fa
+= 3)
656 type
= forceatoms
[fa
];
657 ai
= forceatoms
[fa
+1];
658 aj
= forceatoms
[fa
+2];
661 ki
= pbc_dx_aiuc(pbc
, x
[ai
], x
[aj
], r
);
665 rvec_sub(x
[ai
], x
[aj
], r
);
668 invr
= gmx::invsqrt(r2
);
670 ex
= ip
[type
].orires
.ex
;
671 power
= ip
[type
].orires
.power
;
672 fc
= smooth_fc
*ip
[type
].orires
.kfac
;
673 dev
= od
->otav
[d
] - ip
[type
].orires
.obs
;
676 * there is no real potential when time averaging is applied
678 vtot
+= 0.5*fc
*gmx::square(dev
);
682 /* Calculate the force as the sqrt of tav times instantaneous */
683 devins
= od
->oins
[d
] - ip
[type
].orires
.obs
;
690 dev
= std::sqrt(dev
*devins
);
698 pfac
= fc
*ip
[type
].orires
.c
*invr2
;
699 for (i
= 0; i
< power
; i
++)
703 mvmul(od
->S
[ex
], r
, Sr
);
704 for (i
= 0; i
< DIM
; i
++)
707 -pfac
*dev
*(4*Sr
[i
] - 2*(2+power
)*invr2
*iprod(Sr
, r
)*r
[i
]);
712 ivec_sub(SHIFT_IVEC(g
, ai
), SHIFT_IVEC(g
, aj
), dt
);
716 for (i
= 0; i
< DIM
; i
++)
720 fshift
[ki
][i
] += fij
[i
];
721 fshift
[CENTRAL
][i
] -= fij
[i
];
729 /* Approx. 80*nfa/3 flops */
732 void update_orires_history(t_fcdata
*fcd
, history_t
*hist
)
740 /* Copy the new time averages that have been calculated
741 * in calc_orires_dev.
743 hist
->orire_initf
= od
->exp_min_t_tau
;
744 for (pair
= 0; pair
< od
->nr
; pair
++)
746 for (i
= 0; i
< 5; i
++)
748 hist
->orire_Dtav
[pair
*5+i
] = od
->Dtav
[pair
][i
];