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, 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/legacyheaders/copyrite.h"
44 #include "gromacs/legacyheaders/main.h"
45 #include "gromacs/legacyheaders/network.h"
46 #include "gromacs/legacyheaders/typedefs.h"
47 #include "gromacs/legacyheaders/types/commrec.h"
48 #include "gromacs/legacyheaders/types/fcdata.h"
49 #include "gromacs/legacyheaders/types/ifunc.h"
50 #include "gromacs/legacyheaders/types/mdatom.h"
51 #include "gromacs/linearalgebra/nrjac.h"
52 #include "gromacs/math/do_fit.h"
53 #include "gromacs/math/vec.h"
54 #include "gromacs/pbcutil/ishift.h"
55 #include "gromacs/pbcutil/mshift.h"
56 #include "gromacs/pbcutil/pbc.h"
57 #include "gromacs/topology/mtop_util.h"
58 #include "gromacs/utility/fatalerror.h"
59 #include "gromacs/utility/smalloc.h"
61 void init_orires(FILE *fplog
, const gmx_mtop_t
*mtop
,
64 const t_commrec
*cr
, t_oriresdata
*od
,
67 int i
, j
, d
, ex
, nmol
, *nr_ex
;
70 gmx_mtop_ilistloop_t iloop
;
72 gmx_mtop_atomloop_all_t aloop
;
74 const gmx_multisim_t
*ms
;
76 od
->nr
= gmx_mtop_ftype_count(mtop
, F_ORIRES
);
79 /* Not doing orientation restraints */
85 gmx_fatal(FARGS
, "Orientation restraints do not work with more than one domain (ie. MPI rank).");
87 /* Orientation restraints */
95 od
->fc
= ir
->orires_fc
;
104 iloop
= gmx_mtop_ilistloop_init(mtop
);
105 while (gmx_mtop_ilistloop_next(iloop
, &il
, &nmol
))
107 for (i
= 0; i
< il
[F_ORIRES
].nr
; i
+= 3)
109 ex
= mtop
->ffparams
.iparams
[il
[F_ORIRES
].iatoms
[i
]].orires
.ex
;
113 for (j
= od
->nex
; j
< ex
+1; j
++)
122 snew(od
->S
, od
->nex
);
123 /* When not doing time averaging, the instaneous and time averaged data
124 * are indentical and the pointers can point to the same memory.
126 snew(od
->Dinsl
, od
->nr
);
129 snew(od
->Dins
, od
->nr
);
133 od
->Dins
= od
->Dinsl
;
136 if (ir
->orires_tau
== 0)
144 snew(od
->Dtav
, od
->nr
);
145 od
->edt
= std::exp(-ir
->delta_t
/ir
->orires_tau
);
146 od
->edt_1
= 1.0 - od
->edt
;
148 /* Extend the state with the orires history */
149 state
->flags
|= (1<<estORIRE_INITF
);
150 state
->hist
.orire_initf
= 1;
151 state
->flags
|= (1<<estORIRE_DTAV
);
152 state
->hist
.norire_Dtav
= od
->nr
*5;
153 snew(state
->hist
.orire_Dtav
, state
->hist
.norire_Dtav
);
156 snew(od
->oinsl
, od
->nr
);
159 snew(od
->oins
, od
->nr
);
163 od
->oins
= od
->oinsl
;
165 if (ir
->orires_tau
== 0)
171 snew(od
->otav
, od
->nr
);
173 snew(od
->tmp
, od
->nex
);
174 snew(od
->TMP
, od
->nex
);
175 for (ex
= 0; ex
< od
->nex
; ex
++)
177 snew(od
->TMP
[ex
], 5);
178 for (i
= 0; i
< 5; i
++)
180 snew(od
->TMP
[ex
][i
], 5);
185 for (i
= 0; i
< mtop
->natoms
; i
++)
187 if (ggrpnr(&mtop
->groups
, egcORFIT
, i
) == 0)
192 snew(od
->mref
, od
->nref
);
193 snew(od
->xref
, od
->nref
);
194 snew(od
->xtmp
, od
->nref
);
196 snew(od
->eig
, od
->nex
*12);
198 /* Determine the reference structure on the master node.
199 * Copy it to the other nodes after checking multi compatibility,
200 * so we are sure the subsystems match before copying.
205 aloop
= gmx_mtop_atomloop_all_init(mtop
);
206 while (gmx_mtop_atomloop_all_next(aloop
, &i
, &atom
))
208 if (mtop
->groups
.grpnr
[egcORFIT
] == NULL
||
209 mtop
->groups
.grpnr
[egcORFIT
][i
] == 0)
211 /* Not correct for free-energy with changing masses */
212 od
->mref
[j
] = atom
->m
;
213 if (ms
== NULL
|| MASTERSIM(ms
))
215 copy_rvec(xref
[i
], od
->xref
[j
]);
216 for (d
= 0; d
< DIM
; d
++)
218 com
[d
] += od
->mref
[j
]*xref
[i
][d
];
225 svmul(1.0/mtot
, com
, com
);
226 if (ms
== NULL
|| MASTERSIM(ms
))
228 for (j
= 0; j
< od
->nref
; j
++)
230 rvec_dec(od
->xref
[j
], com
);
234 fprintf(fplog
, "Found %d orientation experiments\n", od
->nex
);
235 for (i
= 0; i
< od
->nex
; i
++)
237 fprintf(fplog
, " experiment %d has %d restraints\n", i
+1, nr_ex
[i
]);
242 fprintf(fplog
, " the fit group consists of %d atoms and has total mass %g\n",
247 fprintf(fplog
, " the orientation restraints are ensemble averaged over %d systems\n", ms
->nsim
);
249 check_multi_int(fplog
, ms
, od
->nr
,
250 "the number of orientation restraints",
252 check_multi_int(fplog
, ms
, od
->nref
,
253 "the number of fit atoms for orientation restraining",
255 check_multi_int(fplog
, ms
, ir
->nsteps
, "nsteps", FALSE
);
256 /* Copy the reference coordinates from the master to the other nodes */
257 gmx_sum_sim(DIM
*od
->nref
, od
->xref
[0], ms
);
260 please_cite(fplog
, "Hess2003");
263 void diagonalize_orires_tensors(t_oriresdata
*od
)
265 int ex
, i
, j
, nrot
, ord
[DIM
], t
;
271 for (i
= 0; i
< DIM
; i
++)
275 snew(od
->eig_diag
, DIM
);
277 for (i
= 0; i
< DIM
; i
++)
283 for (ex
= 0; ex
< od
->nex
; ex
++)
285 /* Rotate the S tensor back to the reference frame */
286 mmul(od
->R
, od
->S
[ex
], TMP
);
287 mtmul(TMP
, od
->R
, S
);
288 for (i
= 0; i
< DIM
; i
++)
290 for (j
= 0; j
< DIM
; j
++)
292 od
->M
[i
][j
] = S
[i
][j
];
296 jacobi(od
->M
, DIM
, od
->eig_diag
, od
->v
, &nrot
);
298 for (i
= 0; i
< DIM
; i
++)
302 for (i
= 0; i
< DIM
; i
++)
304 for (j
= i
+1; j
< DIM
; j
++)
306 if (sqr(od
->eig_diag
[ord
[j
]]) > sqr(od
->eig_diag
[ord
[i
]]))
315 for (i
= 0; i
< DIM
; i
++)
317 od
->eig
[ex
*12 + i
] = od
->eig_diag
[ord
[i
]];
319 for (i
= 0; i
< DIM
; i
++)
321 for (j
= 0; j
< DIM
; j
++)
323 od
->eig
[ex
*12 + 3 + 3*i
+ j
] = od
->v
[j
][ord
[i
]];
329 void print_orires_log(FILE *log
, t_oriresdata
*od
)
334 diagonalize_orires_tensors(od
);
336 for (ex
= 0; ex
< od
->nex
; ex
++)
338 eig
= od
->eig
+ ex
*12;
339 fprintf(log
, " Orientation experiment %d:\n", ex
+1);
340 fprintf(log
, " order parameter: %g\n", eig
[0]);
341 for (i
= 0; i
< DIM
; i
++)
343 fprintf(log
, " eig: %6.3f %6.3f %6.3f %6.3f\n",
344 (eig
[0] != 0) ? eig
[i
]/eig
[0] : eig
[i
],
353 real
calc_orires_dev(const gmx_multisim_t
*ms
,
354 int nfa
, const t_iatom forceatoms
[], const t_iparams ip
[],
355 const t_mdatoms
*md
, const rvec x
[], const t_pbc
*pbc
,
356 t_fcdata
*fcd
, history_t
*hist
)
358 int fa
, d
, i
, j
, type
, ex
, nref
;
359 real edt
, edt_1
, invn
, pfac
, r2
, invr
, corrfac
, weight
, wsv2
, sw
, dev
;
361 rvec5
*Dinsl
, *Dins
, *Dtav
, *rhs
;
364 rvec
*xref
, *xtmp
, com
, r_unrot
, r
;
367 const real two_thr
= 2.0/3.0;
373 /* This means that this is not the master node */
374 gmx_fatal(FARGS
, "Orientation restraints are only supported on the master rank, use fewer ranks");
377 bTAV
= (od
->edt
!= 0);
393 od
->exp_min_t_tau
= hist
->orire_initf
*edt
;
395 /* Correction factor to correct for the lack of history
398 corrfac
= 1.0/(1.0 - od
->exp_min_t_tau
);
417 for (i
= 0; i
< md
->nr
; i
++)
419 if (md
->cORF
[i
] == 0)
421 copy_rvec(x
[i
], xtmp
[j
]);
422 mref
[j
] = md
->massT
[i
];
423 for (d
= 0; d
< DIM
; d
++)
425 com
[d
] += mref
[j
]*xref
[j
][d
];
431 svmul(1.0/mtot
, com
, com
);
432 for (j
= 0; j
< nref
; j
++)
434 rvec_dec(xtmp
[j
], com
);
436 /* Calculate the rotation matrix to rotate x to the reference orientation */
437 calc_fit_R(DIM
, nref
, mref
, xref
, xtmp
, R
);
441 for (fa
= 0; fa
< nfa
; fa
+= 3)
443 type
= forceatoms
[fa
];
446 pbc_dx_aiuc(pbc
, x
[forceatoms
[fa
+1]], x
[forceatoms
[fa
+2]], r_unrot
);
450 rvec_sub(x
[forceatoms
[fa
+1]], x
[forceatoms
[fa
+2]], r_unrot
);
452 mvmul(R
, r_unrot
, r
);
454 invr
= gmx_invsqrt(r2
);
455 /* Calculate the prefactor for the D tensor, this includes the factor 3! */
456 pfac
= ip
[type
].orires
.c
*invr
*invr
*3;
457 for (i
= 0; i
< ip
[type
].orires
.power
; i
++)
461 Dinsl
[d
][0] = pfac
*(2*r
[0]*r
[0] + r
[1]*r
[1] - r2
);
462 Dinsl
[d
][1] = pfac
*(2*r
[0]*r
[1]);
463 Dinsl
[d
][2] = pfac
*(2*r
[0]*r
[2]);
464 Dinsl
[d
][3] = pfac
*(2*r
[1]*r
[1] + r
[0]*r
[0] - r2
);
465 Dinsl
[d
][4] = pfac
*(2*r
[1]*r
[2]);
469 for (i
= 0; i
< 5; i
++)
471 Dins
[d
][i
] = Dinsl
[d
][i
]*invn
;
480 gmx_sum_sim(5*od
->nr
, Dins
[0], ms
);
483 /* Calculate the order tensor S for each experiment via optimization */
484 for (ex
= 0; ex
< od
->nex
; ex
++)
486 for (i
= 0; i
< 5; i
++)
489 for (j
= 0; j
<= i
; j
++)
496 for (fa
= 0; fa
< nfa
; fa
+= 3)
500 /* Here we update Dtav in t_fcdata using the data in history_t.
501 * Thus the results stay correct when this routine
502 * is called multiple times.
504 for (i
= 0; i
< 5; i
++)
506 Dtav
[d
][i
] = edt
*hist
->orire_Dtav
[d
*5+i
] + edt_1
*Dins
[d
][i
];
510 type
= forceatoms
[fa
];
511 ex
= ip
[type
].orires
.ex
;
512 weight
= ip
[type
].orires
.kfac
;
513 /* Calculate the vector rhs and half the matrix T for the 5 equations */
514 for (i
= 0; i
< 5; i
++)
516 rhs
[ex
][i
] += Dtav
[d
][i
]*ip
[type
].orires
.obs
*weight
;
517 for (j
= 0; j
<= i
; j
++)
519 T
[ex
][i
][j
] += Dtav
[d
][i
]*Dtav
[d
][j
]*weight
;
524 /* Now we have all the data we can calculate S */
525 for (ex
= 0; ex
< od
->nex
; ex
++)
527 /* Correct corrfac and copy one half of T to the other half */
528 for (i
= 0; i
< 5; i
++)
530 rhs
[ex
][i
] *= corrfac
;
531 T
[ex
][i
][i
] *= sqr(corrfac
);
532 for (j
= 0; j
< i
; j
++)
534 T
[ex
][i
][j
] *= sqr(corrfac
);
535 T
[ex
][j
][i
] = T
[ex
][i
][j
];
538 m_inv_gen(T
[ex
], 5, T
[ex
]);
539 /* Calculate the orientation tensor S for this experiment */
545 for (i
= 0; i
< 5; i
++)
547 S
[ex
][0][0] += 1.5*T
[ex
][0][i
]*rhs
[ex
][i
];
548 S
[ex
][0][1] += 1.5*T
[ex
][1][i
]*rhs
[ex
][i
];
549 S
[ex
][0][2] += 1.5*T
[ex
][2][i
]*rhs
[ex
][i
];
550 S
[ex
][1][1] += 1.5*T
[ex
][3][i
]*rhs
[ex
][i
];
551 S
[ex
][1][2] += 1.5*T
[ex
][4][i
]*rhs
[ex
][i
];
553 S
[ex
][1][0] = S
[ex
][0][1];
554 S
[ex
][2][0] = S
[ex
][0][2];
555 S
[ex
][2][1] = S
[ex
][1][2];
556 S
[ex
][2][2] = -S
[ex
][0][0] - S
[ex
][1][1];
563 for (fa
= 0; fa
< nfa
; fa
+= 3)
565 type
= forceatoms
[fa
];
566 ex
= ip
[type
].orires
.ex
;
568 od
->otav
[d
] = two_thr
*
569 corrfac
*(S
[ex
][0][0]*Dtav
[d
][0] + S
[ex
][0][1]*Dtav
[d
][1] +
570 S
[ex
][0][2]*Dtav
[d
][2] + S
[ex
][1][1]*Dtav
[d
][3] +
571 S
[ex
][1][2]*Dtav
[d
][4]);
574 od
->oins
[d
] = two_thr
*(S
[ex
][0][0]*Dins
[d
][0] + S
[ex
][0][1]*Dins
[d
][1] +
575 S
[ex
][0][2]*Dins
[d
][2] + S
[ex
][1][1]*Dins
[d
][3] +
576 S
[ex
][1][2]*Dins
[d
][4]);
580 /* When ensemble averaging is used recalculate the local orientation
581 * for output to the energy file.
583 od
->oinsl
[d
] = two_thr
*
584 (S
[ex
][0][0]*Dinsl
[d
][0] + S
[ex
][0][1]*Dinsl
[d
][1] +
585 S
[ex
][0][2]*Dinsl
[d
][2] + S
[ex
][1][1]*Dinsl
[d
][3] +
586 S
[ex
][1][2]*Dinsl
[d
][4]);
589 dev
= od
->otav
[d
] - ip
[type
].orires
.obs
;
591 wsv2
+= ip
[type
].orires
.kfac
*sqr(dev
);
592 sw
+= ip
[type
].orires
.kfac
;
596 od
->rmsdev
= std::sqrt(wsv2
/sw
);
598 /* Rotate the S matrices back, so we get the correct grad(tr(S D)) */
599 for (ex
= 0; ex
< od
->nex
; ex
++)
601 tmmul(R
, S
[ex
], TMP
);
607 /* Approx. 120*nfa/3 flops */
610 real
orires(int nfa
, const t_iatom forceatoms
[], const t_iparams ip
[],
611 const rvec x
[], rvec f
[], rvec fshift
[],
612 const t_pbc
*pbc
, const t_graph
*g
,
613 real gmx_unused lambda
, real gmx_unused
*dvdlambda
,
614 const t_mdatoms gmx_unused
*md
, t_fcdata
*fcd
,
615 int gmx_unused
*global_atom_index
)
618 int fa
, d
, i
, type
, ex
, power
, ki
= CENTRAL
;
620 real r2
, invr
, invr2
, fc
, smooth_fc
, dev
, devins
, pfac
;
623 const t_oriresdata
*od
;
631 bTAV
= (od
->edt
!= 0);
636 /* Smoothly switch on the restraining when time averaging is used */
637 smooth_fc
*= (1.0 - od
->exp_min_t_tau
);
641 for (fa
= 0; fa
< nfa
; fa
+= 3)
643 type
= forceatoms
[fa
];
644 ai
= forceatoms
[fa
+1];
645 aj
= forceatoms
[fa
+2];
648 ki
= pbc_dx_aiuc(pbc
, x
[ai
], x
[aj
], r
);
652 rvec_sub(x
[ai
], x
[aj
], r
);
655 invr
= gmx_invsqrt(r2
);
657 ex
= ip
[type
].orires
.ex
;
658 power
= ip
[type
].orires
.power
;
659 fc
= smooth_fc
*ip
[type
].orires
.kfac
;
660 dev
= od
->otav
[d
] - ip
[type
].orires
.obs
;
663 * there is no real potential when time averaging is applied
665 vtot
+= 0.5*fc
*sqr(dev
);
669 /* Calculate the force as the sqrt of tav times instantaneous */
670 devins
= od
->oins
[d
] - ip
[type
].orires
.obs
;
677 dev
= std::sqrt(dev
*devins
);
685 pfac
= fc
*ip
[type
].orires
.c
*invr2
;
686 for (i
= 0; i
< power
; i
++)
690 mvmul(od
->S
[ex
], r
, Sr
);
691 for (i
= 0; i
< DIM
; i
++)
694 -pfac
*dev
*(4*Sr
[i
] - 2*(2+power
)*invr2
*iprod(Sr
, r
)*r
[i
]);
699 ivec_sub(SHIFT_IVEC(g
, ai
), SHIFT_IVEC(g
, aj
), dt
);
703 for (i
= 0; i
< DIM
; i
++)
707 fshift
[ki
][i
] += fij
[i
];
708 fshift
[CENTRAL
][i
] -= fij
[i
];
716 /* Approx. 80*nfa/3 flops */
719 void update_orires_history(t_fcdata
*fcd
, history_t
*hist
)
727 /* Copy the new time averages that have been calculated
728 * in calc_orires_dev.
730 hist
->orire_initf
= od
->exp_min_t_tau
;
731 for (pair
= 0; pair
< od
->nr
; pair
++)
733 for (i
= 0; i
< 5; i
++)
735 hist
->orire_Dtav
[pair
*5+i
] = od
->Dtav
[pair
][i
];