1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * GROningen Mixture of Alchemy and Childrens' Stories
36 /* This file is completely threadsafe - keep it that way! */
49 #include "gmx_fatal.h"
54 #include "mtop_util.h"
56 void init_disres(FILE *fplog
,const gmx_mtop_t
*mtop
,
57 t_inputrec
*ir
,const t_commrec
*cr
,gmx_bool bPartDecomp
,
58 t_fcdata
*fcd
,t_state
*state
)
60 int fa
,nmol
,i
,npair
,np
;
64 gmx_mtop_ilistloop_t iloop
;
70 if (gmx_mtop_ftype_count(mtop
,F_DISRES
) == 0)
79 fprintf(fplog
,"Initializing the distance restraints\n");
83 if (ir
->eDisre
== edrEnsemble
)
85 gmx_fatal(FARGS
,"Sorry, distance restraints with ensemble averaging over multiple molecules in one system are not functional in this version of GROMACS");
88 dd
->dr_weighting
= ir
->eDisreWeighting
;
89 dd
->dr_fc
= ir
->dr_fc
;
90 if (EI_DYNAMICS(ir
->eI
))
92 dd
->dr_tau
= ir
->dr_tau
;
98 if (dd
->dr_tau
== 0.0)
100 dd
->dr_bMixed
= FALSE
;
105 dd
->dr_bMixed
= ir
->bDisreMixed
;
106 dd
->ETerm
= exp(-(ir
->delta_t
/ir
->dr_tau
));
108 dd
->ETerm1
= 1.0 - dd
->ETerm
;
110 ip
= mtop
->ffparams
.iparams
;
114 iloop
= gmx_mtop_ilistloop_init(mtop
);
115 while (gmx_mtop_ilistloop_next(iloop
,&il
,&nmol
)) {
117 for(fa
=0; fa
<il
[F_DISRES
].nr
; fa
+=3)
120 npair
= mtop
->ffparams
.iparams
[il
[F_DISRES
].iatoms
[fa
]].disres
.npair
;
123 dd
->nres
+= (ir
->eDisre
==edrEnsemble
? 1 : nmol
)*npair
;
124 dd
->npair
+= nmol
*npair
;
130 if (cr
&& PAR(cr
) && !bPartDecomp
)
132 /* Temporary check, will be removed when disre is implemented with DD */
133 const char *notestr
="NOTE: atoms involved in distance restraints should be within the longest cut-off distance, if this is not the case mdrun generates a fatal error, in that case use particle decomposition (mdrun option -pd)";
136 fprintf(stderr
,"\n%s\n\n",notestr
);
138 fprintf(fplog
,"%s\n",notestr
);
140 if (dd
->dr_tau
!= 0 || ir
->eDisre
== edrEnsemble
|| cr
->ms
!= NULL
||
141 dd
->nres
!= dd
->npair
)
143 gmx_fatal(FARGS
,"Time or ensemble averaged or multiple pair distance restraints do not work (yet) with domain decomposition, use particle decomposition (mdrun option -pd)");
145 if (ir
->nstdisreout
!= 0)
149 fprintf(fplog
,"\nWARNING: Can not write distance restraint data to energy file with domain decomposition\n\n");
153 fprintf(stderr
,"\nWARNING: Can not write distance restraint data to energy file with domain decomposition\n");
159 snew(dd
->rt
,dd
->npair
);
161 if (dd
->dr_tau
!= 0.0)
164 /* Set the "history lack" factor to 1 */
165 state
->flags
|= (1<<estDISRE_INITF
);
166 hist
->disre_initf
= 1.0;
167 /* Allocate space for the r^-3 time averages */
168 state
->flags
|= (1<<estDISRE_RM3TAV
);
169 hist
->ndisrepairs
= dd
->npair
;
170 snew(hist
->disre_rm3tav
,hist
->ndisrepairs
);
172 /* Allocate space for a copy of rm3tav,
173 * so we can call do_force without modifying the state.
175 snew(dd
->rm3tav
,dd
->npair
);
177 /* Allocate Rt_6 and Rtav_6 consecutively in memory so they can be
178 * averaged over the processors in one call (in calc_disre_R_6)
180 snew(dd
->Rt_6
,2*dd
->nres
);
181 dd
->Rtav_6
= &(dd
->Rt_6
[dd
->nres
]);
183 ptr
= getenv("GMX_DISRE_ENSEMBLE_SIZE");
184 if (cr
&& cr
->ms
!= NULL
&& ptr
!= NULL
)
188 sscanf(ptr
,"%d",&dd
->nsystems
);
191 fprintf(fplog
,"Found GMX_DISRE_ENSEMBLE_SIZE set to %d systems per ensemble\n",dd
->nsystems
);
193 check_multi_int(fplog
,cr
->ms
,dd
->nsystems
,
194 "the number of systems per ensemble");
195 if (dd
->nsystems
<= 0 || cr
->ms
->nsim
% dd
->nsystems
!= 0)
197 gmx_fatal(FARGS
,"The number of systems %d is not divisible by the number of systems per ensemble %d\n",cr
->ms
->nsim
,dd
->nsystems
);
199 /* Split the inter-master communicator into different ensembles */
200 MPI_Comm_split(cr
->ms
->mpi_comm_masters
,
201 cr
->ms
->sim
/dd
->nsystems
,
203 &dd
->mpi_comm_ensemble
);
206 fprintf(fplog
,"Our ensemble consists of systems:");
207 for(i
=0; i
<dd
->nsystems
; i
++)
210 (cr
->ms
->sim
/dd
->nsystems
)*dd
->nsystems
+i
);
214 snew(dd
->Rtl_6
,dd
->nres
);
220 dd
->Rtl_6
= dd
->Rt_6
;
226 fprintf(fplog
,"There are %d distance restraints involving %d atom pairs\n",dd
->nres
,dd
->npair
);
230 check_multi_int(fplog
,cr
->ms
,fcd
->disres
.nres
,
231 "the number of distance restraints");
233 please_cite(fplog
,"Tropp80a");
234 please_cite(fplog
,"Torda89a");
238 void calc_disres_R_6(const gmx_multisim_t
*ms
,
239 int nfa
,const t_iatom forceatoms
[],const t_iparams ip
[],
240 const rvec x
[],const t_pbc
*pbc
,
241 t_fcdata
*fcd
,history_t
*hist
)
244 int fa
,res
,i
,pair
,ki
,kj
,m
;
247 real
*rt
,*rm3tav
,*Rtl_6
,*Rt_6
,*Rtav_6
;
251 real ETerm
,ETerm1
,cf1
=0,cf2
=0,invn
=0;
255 bTav
= (dd
->dr_tau
!= 0);
266 /* scaling factor to smoothly turn on the restraint forces *
267 * when using time averaging */
268 dd
->exp_min_t_tau
= hist
->disre_initf
*ETerm
;
270 cf1
= dd
->exp_min_t_tau
;
271 cf2
= 1.0/(1.0 - dd
->exp_min_t_tau
);
274 if (dd
->nsystems
> 1)
276 invn
= 1.0/dd
->nsystems
;
279 /* 'loop' over all atom pairs (pair_nr=fa/3) involved in restraints, *
280 * the total number of atoms pairs is nfa/3 */
285 type
= forceatoms
[fa
];
286 npair
= ip
[type
].disres
.npair
;
291 /* Loop over the atom pairs of 'this' restraint */
293 while (fa
< nfa
&& np
< npair
)
296 ai
= forceatoms
[fa
+1];
297 aj
= forceatoms
[fa
+2];
301 pbc_dx_aiuc(pbc
,x
[ai
],x
[aj
],dx
);
305 rvec_sub(x
[ai
],x
[aj
],dx
);
308 rt_1
= gmx_invsqrt(rt2
);
309 rt_3
= rt_1
*rt_1
*rt_1
;
311 rt
[pair
] = sqrt(rt2
);
314 /* Here we update rm3tav in t_fcdata using the data
316 * Thus the results stay correct when this routine
317 * is called multiple times.
319 rm3tav
[pair
] = cf2
*((ETerm
- cf1
)*hist
->disre_rm3tav
[pair
] +
327 Rt_6
[res
] += rt_3
*rt_3
;
328 Rtav_6
[res
] += rm3tav
[pair
]*rm3tav
[pair
];
333 if (dd
->nsystems
> 1)
335 Rtl_6
[res
] = Rt_6
[res
];
344 if (dd
->nsystems
> 1)
346 gmx_sum_comm(2*dd
->nres
,Rt_6
,dd
->mpi_comm_ensemble
);
351 real
ta_disres(int nfa
,const t_iatom forceatoms
[],const t_iparams ip
[],
352 const rvec x
[],rvec f
[],rvec fshift
[],
353 const t_pbc
*pbc
,const t_graph
*g
,
354 real lambda
,real
*dvdlambda
,
355 const t_mdatoms
*md
,t_fcdata
*fcd
,
356 int *global_atom_index
)
358 const real sixth
=1.0/6.0;
359 const real seven_three
=7.0/3.0;
362 int fa
,res
,npair
,p
,pair
,ki
=CENTRAL
,m
;
366 real smooth_fc
,Rt
,Rtav
,rt2
,*Rtl_6
,*Rt_6
,*Rtav_6
;
367 real k0
,f_scal
=0,fmax_scal
,fk_scal
,fij
;
368 real tav_viol
,instant_viol
,mixed_viol
,violtot
,vtot
;
369 real tav_viol_Rtav7
,instant_viol_Rtav7
;
371 gmx_bool bConservative
,bMixed
,bViolation
;
378 dr_weighting
= dd
->dr_weighting
;
379 dr_bMixed
= dd
->dr_bMixed
;
384 tav_viol
=instant_viol
=mixed_viol
=tav_viol_Rtav7
=instant_viol_Rtav7
=0;
386 smooth_fc
= dd
->dr_fc
;
389 /* scaling factor to smoothly turn on the restraint forces *
390 * when using time averaging */
391 smooth_fc
*= (1.0 - dd
->exp_min_t_tau
);
397 /* 'loop' over all atom pairs (pair_nr=fa/3) involved in restraints, *
398 * the total number of atoms pairs is nfa/3 */
403 type
= forceatoms
[fa
];
404 /* Take action depending on restraint, calculate scalar force */
405 npair
= ip
[type
].disres
.npair
;
406 up1
= ip
[type
].disres
.up1
;
407 up2
= ip
[type
].disres
.up2
;
408 low
= ip
[type
].disres
.low
;
409 k0
= smooth_fc
*ip
[type
].disres
.kfac
;
411 /* save some flops when there is only one pair */
412 if (ip
[type
].disres
.type
!= 2)
414 bConservative
= (dr_weighting
== edrwConservative
) && (npair
> 1);
416 Rt
= pow(Rt_6
[res
],-sixth
);
417 Rtav
= pow(Rtav_6
[res
],-sixth
);
421 /* When rtype=2 use instantaneous not ensemble avereged distance */
422 bConservative
= (npair
> 1);
424 Rt
= pow(Rtl_6
[res
],-sixth
);
431 tav_viol
= Rtav
- up1
;
436 tav_viol
= Rtav
- low
;
446 * there is no real potential when time averaging is applied
448 vtot
+= 0.5*k0
*sqr(tav_viol
);
451 printf("vtot is inf: %f\n",vtot
);
455 f_scal
= -k0
*tav_viol
;
456 violtot
+= fabs(tav_viol
);
464 instant_viol
= Rt
- up1
;
475 instant_viol
= Rt
- low
;
488 mixed_viol
= sqrt(tav_viol
*instant_viol
);
489 f_scal
= -k0
*mixed_viol
;
490 violtot
+= mixed_viol
;
497 fmax_scal
= -k0
*(up2
-up1
);
498 /* Correct the force for the number of restraints */
501 f_scal
= max(f_scal
,fmax_scal
);
504 f_scal
*= Rtav
/Rtav_6
[res
];
508 f_scal
/= 2*mixed_viol
;
509 tav_viol_Rtav7
= tav_viol
*Rtav
/Rtav_6
[res
];
510 instant_viol_Rtav7
= instant_viol
*Rt
/Rt_6
[res
];
515 f_scal
/= (real
)npair
;
516 f_scal
= max(f_scal
,fmax_scal
);
519 /* Exert the force ... */
521 /* Loop over the atom pairs of 'this' restraint */
522 for(p
=0; p
<npair
; p
++)
525 ai
= forceatoms
[fa
+1];
526 aj
= forceatoms
[fa
+2];
530 ki
= pbc_dx_aiuc(pbc
,x
[ai
],x
[aj
],dx
);
534 rvec_sub(x
[ai
],x
[aj
],dx
);
538 weight_rt_1
= gmx_invsqrt(rt2
);
544 weight_rt_1
*= pow(dd
->rm3tav
[pair
],seven_three
);
548 weight_rt_1
*= tav_viol_Rtav7
*pow(dd
->rm3tav
[pair
],seven_three
)+
549 instant_viol_Rtav7
*pow(dd
->rt
[pair
],-7);
553 fk_scal
= f_scal
*weight_rt_1
;
557 ivec_sub(SHIFT_IVEC(g
,ai
),SHIFT_IVEC(g
,aj
),dt
);
567 fshift
[ki
][m
] += fij
;
568 fshift
[CENTRAL
][m
] -= fij
;
575 /* No violation so force and potential contributions */
581 dd
->sumviol
= violtot
;
587 void update_disres_history(t_fcdata
*fcd
,history_t
*hist
)
595 /* Copy the new time averages that have been calculated
596 * in calc_disres_R_6.
598 hist
->disre_initf
= dd
->exp_min_t_tau
;
599 for(pair
=0; pair
<dd
->npair
; pair
++)
601 hist
->disre_rm3tav
[pair
] = dd
->rm3tav
[pair
];