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, 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.
41 #include "gromacs/legacyheaders/typedefs.h"
42 #include "gromacs/legacyheaders/types/commrec.h"
43 #include "gromacs/legacyheaders/network.h"
44 #include "gromacs/legacyheaders/orires.h"
45 #include "gromacs/legacyheaders/main.h"
46 #include "gromacs/legacyheaders/copyrite.h"
47 #include "gromacs/pbcutil/pbc.h"
48 #include "gromacs/topology/mtop_util.h"
50 #include "gromacs/linearalgebra/nrjac.h"
51 #include "gromacs/math/do_fit.h"
52 #include "gromacs/math/vec.h"
53 #include "gromacs/pbcutil/ishift.h"
54 #include "gromacs/pbcutil/mshift.h"
55 #include "gromacs/utility/fatalerror.h"
56 #include "gromacs/utility/smalloc.h"
58 void init_orires(FILE *fplog, const gmx_mtop_t *mtop,
61 const t_commrec *cr, t_oriresdata *od,
64 int i, j, d, ex, nmol, *nr_ex;
67 gmx_mtop_ilistloop_t iloop;
69 gmx_mtop_atomloop_all_t aloop;
71 const gmx_multisim_t *ms;
73 od->nr = gmx_mtop_ftype_count(mtop, F_ORIRES);
76 /* Not doing orientation restraints */
82 gmx_fatal(FARGS, "Orientation restraints do not work with more than one domain (ie. MPI rank).");
84 /* Orientation restraints */
92 od->fc = ir->orires_fc;
101 iloop = gmx_mtop_ilistloop_init(mtop);
102 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
104 for (i = 0; i < il[F_ORIRES].nr; i += 3)
106 ex = mtop->ffparams.iparams[il[F_ORIRES].iatoms[i]].orires.ex;
110 for (j = od->nex; j < ex+1; j++)
119 snew(od->S, od->nex);
120 /* When not doing time averaging, the instaneous and time averaged data
121 * are indentical and the pointers can point to the same memory.
123 snew(od->Dinsl, od->nr);
126 snew(od->Dins, od->nr);
130 od->Dins = od->Dinsl;
133 if (ir->orires_tau == 0)
141 snew(od->Dtav, od->nr);
142 od->edt = exp(-ir->delta_t/ir->orires_tau);
143 od->edt_1 = 1.0 - od->edt;
145 /* Extend the state with the orires history */
146 state->flags |= (1<<estORIRE_INITF);
147 state->hist.orire_initf = 1;
148 state->flags |= (1<<estORIRE_DTAV);
149 state->hist.norire_Dtav = od->nr*5;
150 snew(state->hist.orire_Dtav, state->hist.norire_Dtav);
153 snew(od->oinsl, od->nr);
156 snew(od->oins, od->nr);
160 od->oins = od->oinsl;
162 if (ir->orires_tau == 0)
168 snew(od->otav, od->nr);
170 snew(od->tmp, od->nex);
171 snew(od->TMP, od->nex);
172 for (ex = 0; ex < od->nex; ex++)
174 snew(od->TMP[ex], 5);
175 for (i = 0; i < 5; i++)
177 snew(od->TMP[ex][i], 5);
182 for (i = 0; i < mtop->natoms; i++)
184 if (ggrpnr(&mtop->groups, egcORFIT, i) == 0)
189 snew(od->mref, od->nref);
190 snew(od->xref, od->nref);
191 snew(od->xtmp, od->nref);
193 snew(od->eig, od->nex*12);
195 /* Determine the reference structure on the master node.
196 * Copy it to the other nodes after checking multi compatibility,
197 * so we are sure the subsystems match before copying.
202 aloop = gmx_mtop_atomloop_all_init(mtop);
203 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
205 if (mtop->groups.grpnr[egcORFIT] == NULL ||
206 mtop->groups.grpnr[egcORFIT][i] == 0)
208 /* Not correct for free-energy with changing masses */
209 od->mref[j] = atom->m;
210 if (ms == NULL || MASTERSIM(ms))
212 copy_rvec(xref[i], od->xref[j]);
213 for (d = 0; d < DIM; d++)
215 com[d] += od->mref[j]*xref[i][d];
222 svmul(1.0/mtot, com, com);
223 if (ms == NULL || MASTERSIM(ms))
225 for (j = 0; j < od->nref; j++)
227 rvec_dec(od->xref[j], com);
231 fprintf(fplog, "Found %d orientation experiments\n", od->nex);
232 for (i = 0; i < od->nex; i++)
234 fprintf(fplog, " experiment %d has %d restraints\n", i+1, nr_ex[i]);
239 fprintf(fplog, " the fit group consists of %d atoms and has total mass %g\n",
244 fprintf(fplog, " the orientation restraints are ensemble averaged over %d systems\n", ms->nsim);
246 check_multi_int(fplog, ms, od->nr,
247 "the number of orientation restraints",
249 check_multi_int(fplog, ms, od->nref,
250 "the number of fit atoms for orientation restraining",
252 check_multi_int(fplog, ms, ir->nsteps, "nsteps", FALSE);
253 /* Copy the reference coordinates from the master to the other nodes */
254 gmx_sum_sim(DIM*od->nref, od->xref[0], ms);
257 please_cite(fplog, "Hess2003");
260 void diagonalize_orires_tensors(t_oriresdata *od)
262 int ex, i, j, nrot, ord[DIM], t;
268 for (i = 0; i < DIM; i++)
272 snew(od->eig_diag, DIM);
274 for (i = 0; i < DIM; i++)
280 for (ex = 0; ex < od->nex; ex++)
282 /* Rotate the S tensor back to the reference frame */
283 mmul(od->R, od->S[ex], TMP);
284 mtmul(TMP, od->R, S);
285 for (i = 0; i < DIM; i++)
287 for (j = 0; j < DIM; j++)
289 od->M[i][j] = S[i][j];
293 jacobi(od->M, DIM, od->eig_diag, od->v, &nrot);
295 for (i = 0; i < DIM; i++)
299 for (i = 0; i < DIM; i++)
301 for (j = i+1; j < DIM; j++)
303 if (sqr(od->eig_diag[ord[j]]) > sqr(od->eig_diag[ord[i]]))
312 for (i = 0; i < DIM; i++)
314 od->eig[ex*12 + i] = od->eig_diag[ord[i]];
316 for (i = 0; i < DIM; i++)
318 for (j = 0; j < DIM; j++)
320 od->eig[ex*12 + 3 + 3*i + j] = od->v[j][ord[i]];
326 void print_orires_log(FILE *log, t_oriresdata *od)
331 diagonalize_orires_tensors(od);
333 for (ex = 0; ex < od->nex; ex++)
335 eig = od->eig + ex*12;
336 fprintf(log, " Orientation experiment %d:\n", ex+1);
337 fprintf(log, " order parameter: %g\n", eig[0]);
338 for (i = 0; i < DIM; i++)
340 fprintf(log, " eig: %6.3f %6.3f %6.3f %6.3f\n",
341 (eig[0] != 0) ? eig[i]/eig[0] : eig[i],
350 real calc_orires_dev(const gmx_multisim_t *ms,
351 int nfa, const t_iatom forceatoms[], const t_iparams ip[],
352 const t_mdatoms *md, const rvec x[], const t_pbc *pbc,
353 t_fcdata *fcd, history_t *hist)
355 int fa, d, i, j, type, ex, nref;
356 real edt, edt_1, invn, pfac, r2, invr, corrfac, weight, wsv2, sw, dev;
358 rvec5 *Dinsl, *Dins, *Dtav, *rhs;
361 rvec *xref, *xtmp, com, r_unrot, r;
364 const real two_thr = 2.0/3.0;
370 /* This means that this is not the master node */
371 gmx_fatal(FARGS, "Orientation restraints are only supported on the master rank, use fewer ranks");
374 bTAV = (od->edt != 0);
390 od->exp_min_t_tau = hist->orire_initf*edt;
392 /* Correction factor to correct for the lack of history
395 corrfac = 1.0/(1.0 - od->exp_min_t_tau);
414 for (i = 0; i < md->nr; i++)
416 if (md->cORF[i] == 0)
418 copy_rvec(x[i], xtmp[j]);
419 mref[j] = md->massT[i];
420 for (d = 0; d < DIM; d++)
422 com[d] += mref[j]*xref[j][d];
428 svmul(1.0/mtot, com, com);
429 for (j = 0; j < nref; j++)
431 rvec_dec(xtmp[j], com);
433 /* Calculate the rotation matrix to rotate x to the reference orientation */
434 calc_fit_R(DIM, nref, mref, xref, xtmp, R);
438 for (fa = 0; fa < nfa; fa += 3)
440 type = forceatoms[fa];
443 pbc_dx_aiuc(pbc, x[forceatoms[fa+1]], x[forceatoms[fa+2]], r_unrot);
447 rvec_sub(x[forceatoms[fa+1]], x[forceatoms[fa+2]], r_unrot);
449 mvmul(R, r_unrot, r);
451 invr = gmx_invsqrt(r2);
452 /* Calculate the prefactor for the D tensor, this includes the factor 3! */
453 pfac = ip[type].orires.c*invr*invr*3;
454 for (i = 0; i < ip[type].orires.power; i++)
458 Dinsl[d][0] = pfac*(2*r[0]*r[0] + r[1]*r[1] - r2);
459 Dinsl[d][1] = pfac*(2*r[0]*r[1]);
460 Dinsl[d][2] = pfac*(2*r[0]*r[2]);
461 Dinsl[d][3] = pfac*(2*r[1]*r[1] + r[0]*r[0] - r2);
462 Dinsl[d][4] = pfac*(2*r[1]*r[2]);
466 for (i = 0; i < 5; i++)
468 Dins[d][i] = Dinsl[d][i]*invn;
477 gmx_sum_sim(5*od->nr, Dins[0], ms);
480 /* Calculate the order tensor S for each experiment via optimization */
481 for (ex = 0; ex < od->nex; ex++)
483 for (i = 0; i < 5; i++)
486 for (j = 0; j <= i; j++)
493 for (fa = 0; fa < nfa; fa += 3)
497 /* Here we update Dtav in t_fcdata using the data in history_t.
498 * Thus the results stay correct when this routine
499 * is called multiple times.
501 for (i = 0; i < 5; i++)
503 Dtav[d][i] = edt*hist->orire_Dtav[d*5+i] + edt_1*Dins[d][i];
507 type = forceatoms[fa];
508 ex = ip[type].orires.ex;
509 weight = ip[type].orires.kfac;
510 /* Calculate the vector rhs and half the matrix T for the 5 equations */
511 for (i = 0; i < 5; i++)
513 rhs[ex][i] += Dtav[d][i]*ip[type].orires.obs*weight;
514 for (j = 0; j <= i; j++)
516 T[ex][i][j] += Dtav[d][i]*Dtav[d][j]*weight;
521 /* Now we have all the data we can calculate S */
522 for (ex = 0; ex < od->nex; ex++)
524 /* Correct corrfac and copy one half of T to the other half */
525 for (i = 0; i < 5; i++)
527 rhs[ex][i] *= corrfac;
528 T[ex][i][i] *= sqr(corrfac);
529 for (j = 0; j < i; j++)
531 T[ex][i][j] *= sqr(corrfac);
532 T[ex][j][i] = T[ex][i][j];
535 m_inv_gen(T[ex], 5, T[ex]);
536 /* Calculate the orientation tensor S for this experiment */
542 for (i = 0; i < 5; i++)
544 S[ex][0][0] += 1.5*T[ex][0][i]*rhs[ex][i];
545 S[ex][0][1] += 1.5*T[ex][1][i]*rhs[ex][i];
546 S[ex][0][2] += 1.5*T[ex][2][i]*rhs[ex][i];
547 S[ex][1][1] += 1.5*T[ex][3][i]*rhs[ex][i];
548 S[ex][1][2] += 1.5*T[ex][4][i]*rhs[ex][i];
550 S[ex][1][0] = S[ex][0][1];
551 S[ex][2][0] = S[ex][0][2];
552 S[ex][2][1] = S[ex][1][2];
553 S[ex][2][2] = -S[ex][0][0] - S[ex][1][1];
560 for (fa = 0; fa < nfa; fa += 3)
562 type = forceatoms[fa];
563 ex = ip[type].orires.ex;
565 od->otav[d] = two_thr*
566 corrfac*(S[ex][0][0]*Dtav[d][0] + S[ex][0][1]*Dtav[d][1] +
567 S[ex][0][2]*Dtav[d][2] + S[ex][1][1]*Dtav[d][3] +
568 S[ex][1][2]*Dtav[d][4]);
571 od->oins[d] = two_thr*(S[ex][0][0]*Dins[d][0] + S[ex][0][1]*Dins[d][1] +
572 S[ex][0][2]*Dins[d][2] + S[ex][1][1]*Dins[d][3] +
573 S[ex][1][2]*Dins[d][4]);
577 /* When ensemble averaging is used recalculate the local orientation
578 * for output to the energy file.
580 od->oinsl[d] = two_thr*
581 (S[ex][0][0]*Dinsl[d][0] + S[ex][0][1]*Dinsl[d][1] +
582 S[ex][0][2]*Dinsl[d][2] + S[ex][1][1]*Dinsl[d][3] +
583 S[ex][1][2]*Dinsl[d][4]);
586 dev = od->otav[d] - ip[type].orires.obs;
588 wsv2 += ip[type].orires.kfac*sqr(dev);
589 sw += ip[type].orires.kfac;
593 od->rmsdev = sqrt(wsv2/sw);
595 /* Rotate the S matrices back, so we get the correct grad(tr(S D)) */
596 for (ex = 0; ex < od->nex; ex++)
598 tmmul(R, S[ex], TMP);
604 /* Approx. 120*nfa/3 flops */
607 real orires(int nfa, const t_iatom forceatoms[], const t_iparams ip[],
608 const rvec x[], rvec f[], rvec fshift[],
609 const t_pbc *pbc, const t_graph *g,
610 real gmx_unused lambda, real gmx_unused *dvdlambda,
611 const t_mdatoms gmx_unused *md, t_fcdata *fcd,
612 int gmx_unused *global_atom_index)
615 int fa, d, i, type, ex, power, ki = CENTRAL;
617 real r2, invr, invr2, fc, smooth_fc, dev, devins, pfac;
620 const t_oriresdata *od;
628 bTAV = (od->edt != 0);
633 /* Smoothly switch on the restraining when time averaging is used */
634 smooth_fc *= (1.0 - od->exp_min_t_tau);
638 for (fa = 0; fa < nfa; fa += 3)
640 type = forceatoms[fa];
641 ai = forceatoms[fa+1];
642 aj = forceatoms[fa+2];
645 ki = pbc_dx_aiuc(pbc, x[ai], x[aj], r);
649 rvec_sub(x[ai], x[aj], r);
652 invr = gmx_invsqrt(r2);
654 ex = ip[type].orires.ex;
655 power = ip[type].orires.power;
656 fc = smooth_fc*ip[type].orires.kfac;
657 dev = od->otav[d] - ip[type].orires.obs;
660 * there is no real potential when time averaging is applied
662 vtot += 0.5*fc*sqr(dev);
666 /* Calculate the force as the sqrt of tav times instantaneous */
667 devins = od->oins[d] - ip[type].orires.obs;
674 dev = sqrt(dev*devins);
682 pfac = fc*ip[type].orires.c*invr2;
683 for (i = 0; i < power; i++)
687 mvmul(od->S[ex], r, Sr);
688 for (i = 0; i < DIM; i++)
691 -pfac*dev*(4*Sr[i] - 2*(2+power)*invr2*iprod(Sr, r)*r[i]);
696 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
700 for (i = 0; i < DIM; i++)
704 fshift[ki][i] += fij[i];
705 fshift[CENTRAL][i] -= fij[i];
713 /* Approx. 80*nfa/3 flops */
716 void update_orires_history(t_fcdata *fcd, history_t *hist)
724 /* Copy the new time averages that have been calculated
725 * in calc_orires_dev.
727 hist->orire_initf = od->exp_min_t_tau;
728 for (pair = 0; pair < od->nr; pair++)
730 for (i = 0; i < 5; i++)
732 hist->orire_Dtav[pair*5+i] = od->Dtav[pair][i];