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 * GROwing Monsters And Cloning Shrimps
54 #include "mtop_util.h"
57 static bool bEInd
[egNR
] = { TRUE
, TRUE
, FALSE
, FALSE
, FALSE
, FALSE
, FALSE
};
59 static bool bEner
[F_NRE
];
61 static const char *conrmsd_nm
[] = { "Constr. rmsd", "Constr.2 rmsd" };
63 static const char *boxs_nm
[] = { "Box-X", "Box-Y", "Box-Z" };
65 static const char *tricl_boxs_nm
[] = { "Box-XX", "Box-YX", "Box-YY",
66 "Box-ZX", "Box-ZY", "Box-ZZ" };
68 static const char *vol_nm
[] = { "Volume" };
70 static const char *dens_nm
[] = {"Density" };
72 static const char *pv_nm
[] = {"pV" };
74 static const char *boxvel_nm
[] = {
75 "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
76 "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
79 #define NBOXS asize(boxs_nm)
80 #define NTRICLBOXS asize(tricl_boxs_nm)
82 static bool bConstr
,bConstrVir
,bTricl
,bDynBox
;
83 static int f_nre
=0,epc
,etc
,nCrmsd
;
85 t_mdebin
*init_mdebin(int fp_ene
,
86 const gmx_mtop_t
*mtop
,
89 const char *ener_nm
[F_NRE
];
90 static const char *vir_nm
[] = {
91 "Vir-XX", "Vir-XY", "Vir-XZ",
92 "Vir-YX", "Vir-YY", "Vir-YZ",
93 "Vir-ZX", "Vir-ZY", "Vir-ZZ"
95 static const char *sv_nm
[] = {
96 "ShakeVir-XX", "ShakeVir-XY", "ShakeVir-XZ",
97 "ShakeVir-YX", "ShakeVir-YY", "ShakeVir-YZ",
98 "ShakeVir-ZX", "ShakeVir-ZY", "ShakeVir-ZZ"
100 static const char *fv_nm
[] = {
101 "ForceVir-XX", "ForceVir-XY", "ForceVir-XZ",
102 "ForceVir-YX", "ForceVir-YY", "ForceVir-YZ",
103 "ForceVir-ZX", "ForceVir-ZY", "ForceVir-ZZ"
105 static const char *pres_nm
[] = {
106 "Pres-XX","Pres-XY","Pres-XZ",
107 "Pres-YX","Pres-YY","Pres-YZ",
108 "Pres-ZX","Pres-ZY","Pres-ZZ"
110 static const char *surft_nm
[] = {
113 static const char *mu_nm
[] = {
114 "Mu-X", "Mu-Y", "Mu-Z"
116 static const char *vcos_nm
[] = {
119 static const char *visc_nm
[] = {
123 const gmx_groups_t
*groups
;
127 int i
,j
,ni
,nj
,n
,k
,kk
,ncon
,nset
;
130 groups
= &mtop
->groups
;
132 bBHAM
= (mtop
->ffparams
.functype
[0] == F_BHAM
);
133 b14
= (gmx_mtop_ftype_count(mtop
,F_LJ14
) > 0 ||
134 gmx_mtop_ftype_count(mtop
,F_LJC14_Q
) > 0);
136 ncon
= gmx_mtop_ftype_count(mtop
,F_CONSTR
);
137 nset
= gmx_mtop_ftype_count(mtop
,F_SETTLE
);
138 bConstr
= (ncon
> 0 || nset
> 0);
141 if (ncon
> 0 && ir
->eConstrAlg
== econtLINCS
) {
147 bConstrVir
= (getenv("GMX_CONSTRAINTVIR") != NULL
);
152 for(i
=0; i
<F_NRE
; i
++) {
156 else if (i
== F_BHAM
)
159 bEner
[i
] = ir
->bQMMM
;
160 else if (i
== F_COUL_LR
)
161 bEner
[i
] = (ir
->rcoulomb
> ir
->rlist
);
162 else if (i
== F_LJ_LR
)
163 bEner
[i
] = (!bBHAM
&& ir
->rvdw
> ir
->rlist
);
164 else if (i
== F_BHAM_LR
)
165 bEner
[i
] = (bBHAM
&& ir
->rvdw
> ir
->rlist
);
166 else if (i
== F_RF_EXCL
)
167 bEner
[i
] = (EEL_RF(ir
->coulombtype
) && ir
->coulombtype
!= eelRF_NEC
);
168 else if (i
== F_COUL_RECIP
)
169 bEner
[i
] = EEL_FULL(ir
->coulombtype
);
170 else if (i
== F_LJ14
)
172 else if (i
== F_COUL14
)
174 else if (i
== F_LJC14_Q
|| i
== F_LJC_PAIRS_NB
)
176 else if ((i
== F_DVDL
) || (i
== F_DKDL
))
177 bEner
[i
] = (ir
->efep
!= efepNO
);
178 else if (i
== F_DHDL_CON
)
179 bEner
[i
] = (ir
->efep
!= efepNO
&& bConstr
);
180 else if ((interaction_function
[i
].flags
& IF_VSITE
) ||
181 (i
== F_CONSTR
) || (i
== F_CONSTRNC
) || (i
== F_SETTLE
))
183 else if ((i
== F_COUL_SR
) || (i
== F_EPOT
) || (i
== F_PRES
) || (i
==F_EQM
))
185 else if ((i
== F_ETOT
) || (i
== F_EKIN
) || (i
== F_TEMP
))
186 bEner
[i
] = EI_DYNAMICS(ir
->eI
);
187 else if (i
== F_DISPCORR
|| i
== F_PDISPCORR
)
188 bEner
[i
] = (ir
->eDispCorr
!= edispcNO
);
189 else if (i
== F_DISRESVIOL
)
190 bEner
[i
] = (gmx_mtop_ftype_count(mtop
,F_DISRES
) > 0);
191 else if (i
== F_ORIRESDEV
)
192 bEner
[i
] = (gmx_mtop_ftype_count(mtop
,F_ORIRES
) > 0);
193 else if (i
== F_CONNBONDS
)
195 else if (i
== F_COM_PULL
)
196 bEner
[i
] = (ir
->ePull
== epullUMBRELLA
|| ir
->ePull
== epullCONST_F
|| ir
->bRot
);
197 else if (i
== F_ECONSERVED
)
198 bEner
[i
] = ((ir
->etc
== etcNOSEHOOVER
|| ir
->etc
== etcVRESCALE
) &&
201 bEner
[i
] = (gmx_mtop_ftype_count(mtop
,i
) > 0);
204 for(i
=0; i
<F_NRE
; i
++)
208 /* FIXME: The constness should not be cast away */
209 ener_nm
[f_nre
]=(char *)interaction_function
[i
].longname
;
215 bTricl
= TRICLINIC(ir
->compress
) || TRICLINIC(ir
->deform
);
216 bDynBox
= DYNAMIC_BOX(*ir
);
219 /* Energy monitoring */
221 md
->ebin
= mk_ebin();
222 /* Pass NULL for unit to let get_ebin_space determine the units
223 * for interaction_function[i].longname
225 md
->ie
= get_ebin_space(md
->ebin
,f_nre
,ener_nm
,NULL
);
228 /* This should be called directly after the call for md->ie,
229 * such that md->iconrmsd follows directly in the list.
231 md
->iconrmsd
= get_ebin_space(md
->ebin
,nCrmsd
,conrmsd_nm
,"");
235 md
->ib
= get_ebin_space(md
->ebin
, bTricl
? NTRICLBOXS
:
236 NBOXS
, bTricl
? tricl_boxs_nm
: boxs_nm
,
238 md
->ivol
= get_ebin_space(md
->ebin
, 1, vol_nm
, unit_volume
);
239 md
->idens
= get_ebin_space(md
->ebin
, 1, dens_nm
, unit_density_SI
);
240 md
->ipv
= get_ebin_space(md
->ebin
, 1, pv_nm
, unit_energy
);
244 md
->isvir
= get_ebin_space(md
->ebin
,asize(sv_nm
),sv_nm
,unit_energy
);
245 md
->ifvir
= get_ebin_space(md
->ebin
,asize(fv_nm
),fv_nm
,unit_energy
);
247 md
->ivir
= get_ebin_space(md
->ebin
,asize(vir_nm
),vir_nm
,unit_energy
);
248 md
->ipres
= get_ebin_space(md
->ebin
,asize(pres_nm
),pres_nm
,unit_pres_bar
);
249 md
->isurft
= get_ebin_space(md
->ebin
,asize(surft_nm
),surft_nm
,
251 if (epc
== epcPARRINELLORAHMAN
)
253 md
->ipc
= get_ebin_space(md
->ebin
,bTricl
? 6 : 3,boxvel_nm
,unit_vel
);
255 md
->imu
= get_ebin_space(md
->ebin
,asize(mu_nm
),mu_nm
,unit_dipole_D
);
256 if (ir
->cos_accel
!= 0)
258 md
->ivcos
= get_ebin_space(md
->ebin
,asize(vcos_nm
),vcos_nm
,unit_vel
);
259 md
->ivisc
= get_ebin_space(md
->ebin
,asize(visc_nm
),visc_nm
,
262 if (ir
->rcoulomb
> ir
->rlist
)
264 bEInd
[egCOULLR
] = TRUE
;
268 if (ir
->rvdw
> ir
->rlist
)
270 bEInd
[egLJLR
] = TRUE
;
275 bEInd
[egLJSR
] = FALSE
;
276 bEInd
[egBHAMSR
] = TRUE
;
277 if (ir
->rvdw
> ir
->rlist
)
279 bEInd
[egBHAMLR
] = TRUE
;
284 bEInd
[egLJ14
] = TRUE
;
285 bEInd
[egCOUL14
] = TRUE
;
288 for(i
=0; (i
<egNR
); i
++)
296 n
=groups
->grps
[egcENER
].nr
;
299 snew(md
->igrp
,md
->nE
);
304 for(k
=0; (k
<md
->nEc
); k
++)
308 for(i
=0; (i
<groups
->grps
[egcENER
].nr
); i
++)
310 ni
=groups
->grps
[egcENER
].nm_ind
[i
];
311 for(j
=i
; (j
<groups
->grps
[egcENER
].nr
); j
++)
313 nj
=groups
->grps
[egcENER
].nm_ind
[j
];
314 for(k
=kk
=0; (k
<egNR
); k
++)
318 sprintf(gnm
[kk
],"%s:%s-%s",egrp_nm
[k
],
319 *(groups
->grpname
[ni
]),*(groups
->grpname
[nj
]));
323 md
->igrp
[n
]=get_ebin_space(md
->ebin
,md
->nEc
,
324 (const char **)gnm
,unit_energy
);
328 for(k
=0; (k
<md
->nEc
); k
++)
336 gmx_incons("Number of energy terms wrong");
340 md
->nTC
=groups
->grps
[egcTC
].nr
;
341 snew(md
->grpnms
,md
->nTC
);
342 snew(md
->tmp_r
,md
->nTC
);
343 snew(md
->tmp_v
,md
->nTC
);
345 for(i
=0; (i
<md
->nTC
); i
++)
347 ni
=groups
->grps
[egcTC
].nm_ind
[i
];
348 sprintf(buf
,"T-%s",*(groups
->grpname
[ni
]));
349 grpnms
[i
]=strdup(buf
);
351 md
->itemp
=get_ebin_space(md
->ebin
,md
->nTC
,(const char **)grpnms
,
353 if (etc
== etcNOSEHOOVER
)
355 for(i
=0; (i
<md
->nTC
); i
++)
357 ni
=groups
->grps
[egcTC
].nm_ind
[i
];
358 sprintf(buf
,"Xi-%s",*(groups
->grpname
[ni
]));
359 grpnms
[i
]=strdup(buf
);
361 md
->itc
=get_ebin_space(md
->ebin
,md
->nTC
,(const char **)grpnms
,
364 else if (etc
== etcBERENDSEN
|| etc
== etcYES
|| etc
== etcVRESCALE
)
366 for(i
=0; (i
<md
->nTC
); i
++)
368 ni
=groups
->grps
[egcTC
].nm_ind
[i
];
369 sprintf(buf
,"Lamb-%s",*(groups
->grpname
[ni
]));
370 grpnms
[i
]=strdup(buf
);
372 md
->itc
=get_ebin_space(md
->ebin
,md
->nTC
,(const char **)grpnms
,"");
376 md
->nU
=groups
->grps
[egcACC
].nr
;
379 snew(grpnms
,3*md
->nU
);
380 for(i
=0; (i
<md
->nU
); i
++)
382 ni
=groups
->grps
[egcACC
].nm_ind
[i
];
383 sprintf(buf
,"Ux-%s",*(groups
->grpname
[ni
]));
384 grpnms
[3*i
+XX
]=strdup(buf
);
385 sprintf(buf
,"Uy-%s",*(groups
->grpname
[ni
]));
386 grpnms
[3*i
+YY
]=strdup(buf
);
387 sprintf(buf
,"Uz-%s",*(groups
->grpname
[ni
]));
388 grpnms
[3*i
+ZZ
]=strdup(buf
);
390 md
->iu
=get_ebin_space(md
->ebin
,3*md
->nU
,(const char **)grpnms
,unit_vel
);
396 do_enxnms(fp_ene
,&md
->ebin
->nener
,&md
->ebin
->enm
);
402 FILE *open_dhdl(const char *filename
,t_inputrec
*ir
)
405 const char *dhdl
="dH/d\\8l\\4",*deltag
="\\8D\\4H",*lambda
="\\8l\\4";
406 char title
[STRLEN
],label_x
[STRLEN
],label_y
[STRLEN
];
408 char **setname
,buf
[STRLEN
];
410 sprintf(label_x
,"%s (%s)","Time",unit_time
);
411 if (ir
->n_flambda
== 0)
413 sprintf(title
,"%s",dhdl
);
414 sprintf(label_y
,"%s (%s %s)",
415 dhdl
,unit_energy
,"[\\8l\\4]\\S-1\\N");
419 sprintf(title
,"%s, %s",dhdl
,deltag
);
420 sprintf(label_y
,"(%s)",unit_energy
);
422 fp
= xvgropen(filename
,title
,label_x
,label_y
);
424 if (ir
->n_flambda
> 0)
426 /* g_bar has to determine the lambda values used in this simulation
427 * from this xvg legend.
429 nsets
= 1 + ir
->n_flambda
;
431 sprintf(buf
,"%s %s %g",dhdl
,lambda
,ir
->init_lambda
);
432 setname
[0] = strdup(buf
);
433 for(s
=1; s
<nsets
; s
++)
435 sprintf(buf
,"%s %s %g",deltag
,lambda
,ir
->flambda
[s
-1]);
436 setname
[s
] = strdup(buf
);
438 xvgr_legend(fp
,nsets
,setname
);
440 for(s
=0; s
<nsets
; s
++)
450 static void copy_energy(real e
[],real ecpy
[])
454 for(i
=j
=0; (i
<F_NRE
); i
++)
458 gmx_incons("Number of energy terms wrong");
461 void upd_mdebin(t_mdebin
*md
,FILE *fp_dhdl
,
465 gmx_enerdata_t
*enerd
,
472 gmx_ekindata_t
*ekind
,
476 int i
,j
,k
,kk
,m
,n
,gid
;
477 real crmsd
[2],tmp6
[6];
478 real bs
[NTRICLBOXS
],vol
,dens
,pv
;
483 /* Do NOT use the box in the state variable, but the separate box provided
484 * as an argument. This is because we sometimes need to write the box from
485 * the last timestep to match the trajectory frames.
487 copy_energy(enerd
->term
,ecopy
);
488 add_ebin(md
->ebin
,md
->ie
,f_nre
,ecopy
,bSum
);
491 crmsd
[0] = constr_rmsd(constr
,FALSE
);
494 crmsd
[1] = constr_rmsd(constr
,TRUE
);
496 add_ebin(md
->ebin
,md
->iconrmsd
,nCrmsd
,crmsd
,FALSE
);
515 vol
= box
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
516 dens
= (tmass
*AMU
)/(vol
*NANO
*NANO
*NANO
);
518 /* This is pV (in kJ/mol) */
519 pv
= vol
*enerd
->term
[F_PRES
]/PRESFAC
;
521 add_ebin(md
->ebin
,md
->ib
,NBOXS
,bs
,bSum
);
522 add_ebin(md
->ebin
,md
->ivol
,1 ,&vol
,bSum
);
523 add_ebin(md
->ebin
,md
->idens
,1 ,&dens
,bSum
);
524 add_ebin(md
->ebin
,md
->ipv
,1 ,&pv
,bSum
);
528 add_ebin(md
->ebin
,md
->isvir
,9,svir
[0],bSum
);
529 add_ebin(md
->ebin
,md
->ifvir
,9,fvir
[0],bSum
);
531 add_ebin(md
->ebin
,md
->ivir
,9,vir
[0],bSum
);
532 add_ebin(md
->ebin
,md
->ipres
,9,pres
[0],bSum
);
533 tmp
= (pres
[ZZ
][ZZ
]-(pres
[XX
][XX
]+pres
[YY
][YY
])*0.5)*box
[ZZ
][ZZ
];
534 add_ebin(md
->ebin
,md
->isurft
,1,&tmp
,bSum
);
535 if (epc
== epcPARRINELLORAHMAN
)
537 tmp6
[0] = state
->boxv
[XX
][XX
];
538 tmp6
[1] = state
->boxv
[YY
][YY
];
539 tmp6
[2] = state
->boxv
[ZZ
][ZZ
];
540 tmp6
[3] = state
->boxv
[YY
][XX
];
541 tmp6
[4] = state
->boxv
[ZZ
][XX
];
542 tmp6
[5] = state
->boxv
[ZZ
][YY
];
543 add_ebin(md
->ebin
,md
->ipc
,bTricl
? 6 : 3,tmp6
,bSum
);
545 add_ebin(md
->ebin
,md
->imu
,3,mu_tot
,bSum
);
546 if (ekind
&& ekind
->cosacc
.cos_accel
!= 0)
548 vol
= box
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
549 dens
= (tmass
*AMU
)/(vol
*NANO
*NANO
*NANO
);
550 add_ebin(md
->ebin
,md
->ivcos
,1,&(ekind
->cosacc
.vcos
),bSum
);
551 /* 1/viscosity, unit 1/(kg m^-1 s^-1) */
552 tmp
= 1/(ekind
->cosacc
.cos_accel
/(ekind
->cosacc
.vcos
*PICO
)
553 *vol
*sqr(box
[ZZ
][ZZ
]*NANO
/(2*M_PI
)));
554 add_ebin(md
->ebin
,md
->ivisc
,1,&tmp
,bSum
);
559 for(i
=0; (i
<md
->nEg
); i
++)
561 for(j
=i
; (j
<md
->nEg
); j
++)
563 gid
=GID(i
,j
,md
->nEg
);
564 for(k
=kk
=0; (k
<egNR
); k
++)
568 eee
[kk
++] = enerd
->grpp
.ener
[k
][gid
];
571 add_ebin(md
->ebin
,md
->igrp
[n
],md
->nEc
,eee
,bSum
);
579 for(i
=0; (i
<md
->nTC
); i
++)
581 md
->tmp_r
[i
] = ekind
->tcstat
[i
].T
;
583 add_ebin(md
->ebin
,md
->itemp
,md
->nTC
,md
->tmp_r
,bSum
);
584 if (etc
== etcNOSEHOOVER
)
586 for(i
=0; (i
<md
->nTC
); i
++)
588 md
->tmp_r
[i
] = state
->nosehoover_xi
[i
];
590 add_ebin(md
->ebin
,md
->itc
,md
->nTC
,md
->tmp_r
,bSum
);
592 else if (etc
== etcBERENDSEN
|| etc
== etcYES
|| etc
== etcVRESCALE
)
594 for(i
=0; (i
<md
->nTC
); i
++)
596 md
->tmp_r
[i
] = ekind
->tcstat
[i
].lambda
;
598 add_ebin(md
->ebin
,md
->itc
,md
->nTC
,md
->tmp_r
,bSum
);
602 if (ekind
&& md
->nU
> 1)
604 for(i
=0; (i
<md
->nU
); i
++)
606 copy_rvec(ekind
->grpstat
[i
].u
,md
->tmp_v
[i
]);
608 add_ebin(md
->ebin
,md
->iu
,3*md
->nU
,md
->tmp_v
[0],bSum
);
611 ebin_increase_count(md
->ebin
,bSum
);
615 fprintf(fp_dhdl
,"%.4f %g",
617 enerd
->term
[F_DVDL
]+enerd
->term
[F_DKDL
]+enerd
->term
[F_DHDL_CON
]);
618 for(i
=1; i
<enerd
->n_lambda
; i
++)
620 fprintf(fp_dhdl
," %g",
621 enerd
->enerpart_lambda
[i
]-enerd
->enerpart_lambda
[0]);
623 fprintf(fp_dhdl
,"\n");
627 void upd_mdebin_step(t_mdebin
*md
)
629 ebin_increase_count(md
->ebin
,FALSE
);
632 static void npr(FILE *log
,int n
,char c
)
634 for(; (n
>0); n
--) fprintf(log
,"%c",c
);
637 static void pprint(FILE *log
,const char *s
,t_mdebin
*md
)
641 char buf1
[22],buf2
[22];
644 fprintf(log
,"\t<====== ");
646 fprintf(log
," ==>\n");
647 fprintf(log
,"\t<==== %s ====>\n",s
);
648 fprintf(log
,"\t<== ");
650 fprintf(log
," ======>\n\n");
652 fprintf(log
,"\tStatistics over %s steps using %s frames\n",
653 gmx_step_str(md
->ebin
->nsteps_sim
,buf1
),
654 gmx_step_str(md
->ebin
->nsum_sim
,buf2
));
658 void print_ebin_header(FILE *log
,gmx_step_t steps
,double time
,real lamb
)
662 fprintf(log
," %12s %12s %12s\n"
663 " %12s %12.5f %12.5f\n\n",
664 "Step","Time","Lambda",gmx_step_str(steps
,buf
),time
,lamb
);
667 void print_ebin(int fp_ene
,bool bEne
,bool bDR
,bool bOR
,
669 gmx_step_t step
,double time
,
670 int mode
,bool bCompact
,
671 t_mdebin
*md
,t_fcdata
*fcd
,
672 gmx_groups_t
*groups
,t_grpopts
*opts
)
674 static char **grpnms
=NULL
;
676 int i
,j
,n
,ni
,nj
,ndr
,nor
;
686 fr
.nsteps
= md
->ebin
->nsteps
;
687 fr
.nsum
= md
->ebin
->nsum
;
688 fr
.nre
= (bEne
) ? md
->ebin
->nener
: 0;
689 fr
.ener
= md
->ebin
->e
;
690 fr
.ndisre
= bDR
? fcd
->disres
.npair
: 0;
691 fr
.disre_rm3tav
= fcd
->disres
.rm3tav
;
692 fr
.disre_rt
= fcd
->disres
.rt
;
693 /* Optional additional blocks */
694 for(i
=0; i
<enxNR
; i
++)
698 if (fcd
->orires
.nr
> 0 && bOR
)
700 diagonalize_orires_tensors(&(fcd
->orires
));
701 nr
[enxOR
] = fcd
->orires
.nr
;
702 block
[enxOR
] = fcd
->orires
.otav
;
703 nr
[enxORI
] = (fcd
->orires
.oinsl
!= fcd
->orires
.otav
) ?
705 block
[enxORI
] = fcd
->orires
.oinsl
;
706 nr
[enxORT
] = fcd
->orires
.nex
*12;
707 block
[enxORT
] = fcd
->orires
.eig
;
710 for(i
=0; i
<enxNR
; i
++)
719 if (fr
.nre
|| fr
.ndisre
|| fr
.nr
[enxOR
] || fr
.nr
[enxORI
])
724 /* We have stored the sums, so reset the sum history */
725 reset_ebin_sums(md
->ebin
);
732 pprint(log
,"A V E R A G E S",md
);
738 pprint(log
,"R M S - F L U C T U A T I O N S",md
);
742 gmx_fatal(FARGS
,"Invalid print mode (%d)",mode
);
747 for(i
=0;i
<opts
->ngtc
;i
++)
749 if(opts
->annealing
[i
]!=eannNO
)
751 fprintf(log
,"Current ref_t for group %s: %8.1f\n",
752 *(groups
->grpname
[groups
->grps
[egcTC
].nm_ind
[i
]]),opts
->ref_t
[i
]);
755 if (mode
==eprNORMAL
&& fcd
->orires
.nr
>0)
757 print_orires_log(log
,&(fcd
->orires
));
759 fprintf(log
," Energies (%s)\n",unit_energy
);
760 pr_ebin(log
,md
->ebin
,md
->ie
,f_nre
+nCrmsd
,5,mode
,TRUE
);
767 pr_ebin(log
,md
->ebin
,md
->ib
, bTricl
? NTRICLBOXS
: NBOXS
,5,mode
,TRUE
);
772 fprintf(log
," Constraint Virial (%s)\n",unit_energy
);
773 pr_ebin(log
,md
->ebin
,md
->isvir
,9,3,mode
,FALSE
);
775 fprintf(log
," Force Virial (%s)\n",unit_energy
);
776 pr_ebin(log
,md
->ebin
,md
->ifvir
,9,3,mode
,FALSE
);
779 fprintf(log
," Total Virial (%s)\n",unit_energy
);
780 pr_ebin(log
,md
->ebin
,md
->ivir
,9,3,mode
,FALSE
);
782 fprintf(log
," Pressure (%s)\n",unit_pres_bar
);
783 pr_ebin(log
,md
->ebin
,md
->ipres
,9,3,mode
,FALSE
);
785 fprintf(log
," Total Dipole (%s)\n",unit_dipole_D
);
786 pr_ebin(log
,md
->ebin
,md
->imu
,3,3,mode
,FALSE
);
795 for(i
=0; (i
<md
->nEg
); i
++)
797 ni
=groups
->grps
[egcENER
].nm_ind
[i
];
798 for(j
=i
; (j
<md
->nEg
); j
++)
800 nj
=groups
->grps
[egcENER
].nm_ind
[j
];
801 sprintf(buf
,"%s-%s",*(groups
->grpname
[ni
]),*(groups
->grpname
[nj
]));
802 grpnms
[n
++]=strdup(buf
);
806 sprintf(buf
,"Epot (%s)",unit_energy
);
807 fprintf(log
,"%15s ",buf
);
808 for(i
=0; (i
<egNR
); i
++)
812 fprintf(log
,"%12s ",egrp_nm
[i
]);
816 for(i
=0; (i
<md
->nE
); i
++)
818 fprintf(log
,"%15s",grpnms
[i
]);
819 pr_ebin(log
,md
->ebin
,md
->igrp
[i
],md
->nEc
,md
->nEc
,mode
,FALSE
);
825 pr_ebin(log
,md
->ebin
,md
->itemp
,md
->nTC
,4,mode
,TRUE
);
830 fprintf(log
,"%15s %12s %12s %12s\n",
831 "Group","Ux","Uy","Uz");
832 for(i
=0; (i
<md
->nU
); i
++)
834 ni
=groups
->grps
[egcACC
].nm_ind
[i
];
835 fprintf(log
,"%15s",*groups
->grpname
[ni
]);
836 pr_ebin(log
,md
->ebin
,md
->iu
+3*i
,3,3,mode
,FALSE
);
845 init_energyhistory(energyhistory_t
* enerhist
)
849 enerhist
->ener_ave
= NULL
;
850 enerhist
->ener_sum
= NULL
;
851 enerhist
->ener_sum_sim
= NULL
;
853 enerhist
->nsteps
= 0;
855 enerhist
->nsteps_sim
= 0;
856 enerhist
->nsum_sim
= 0;
860 update_energyhistory(energyhistory_t
* enerhist
,t_mdebin
* mdebin
)
864 enerhist
->nsteps
= mdebin
->ebin
->nsteps
;
865 enerhist
->nsum
= mdebin
->ebin
->nsum
;
866 enerhist
->nsteps_sim
= mdebin
->ebin
->nsteps_sim
;
867 enerhist
->nsum_sim
= mdebin
->ebin
->nsum_sim
;
868 enerhist
->nener
= mdebin
->ebin
->nener
;
870 if (mdebin
->ebin
->nsum
> 0)
872 /* Check if we need to allocate first */
873 if(enerhist
->ener_ave
== NULL
)
875 snew(enerhist
->ener_ave
,enerhist
->nener
);
876 snew(enerhist
->ener_sum
,enerhist
->nener
);
879 for(i
=0;i
<enerhist
->nener
;i
++)
881 enerhist
->ener_ave
[i
] = mdebin
->ebin
->e
[i
].eav
;
882 enerhist
->ener_sum
[i
] = mdebin
->ebin
->e
[i
].esum
;
886 if (mdebin
->ebin
->nsum_sim
> 0)
888 /* Check if we need to allocate first */
889 if(enerhist
->ener_sum_sim
== NULL
)
891 snew(enerhist
->ener_sum_sim
,enerhist
->nener
);
894 for(i
=0;i
<enerhist
->nener
;i
++)
896 enerhist
->ener_sum_sim
[i
] = mdebin
->ebin
->e_sim
[i
].esum
;
902 restore_energyhistory_from_state(t_mdebin
* mdebin
,energyhistory_t
* enerhist
)
906 if ((enerhist
->nsum
> 0 || enerhist
->nsum_sim
> 0) &&
907 mdebin
->ebin
->nener
!= enerhist
->nener
)
909 gmx_fatal(FARGS
,"Mismatch between number of energies in run input (%d) and checkpoint file (%d).",
910 mdebin
->ebin
->nener
,enerhist
->nener
);
913 mdebin
->ebin
->nsteps
= enerhist
->nsteps
;
914 mdebin
->ebin
->nsum
= enerhist
->nsum
;
915 mdebin
->ebin
->nsteps_sim
= enerhist
->nsteps_sim
;
916 mdebin
->ebin
->nsum_sim
= enerhist
->nsum_sim
;
918 for(i
=0; i
<mdebin
->ebin
->nener
; i
++)
920 mdebin
->ebin
->e
[i
].eav
=
921 (enerhist
->nsum
> 0 ? enerhist
->ener_ave
[i
] : 0);
922 mdebin
->ebin
->e
[i
].esum
=
923 (enerhist
->nsum
> 0 ? enerhist
->ener_sum
[i
] : 0);
924 mdebin
->ebin
->e_sim
[i
].esum
=
925 (enerhist
->nsum_sim
> 0 ? enerhist
->ener_sum_sim
[i
] : 0);