?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Improved version for builds from modified trees.
[gromacs/qmmm-gamess-us.git] / src / mdlib / mdebin.c
blobeac5dff2eb1a71779b2fac283af088cc36bc8232
1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
2 *
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 3.2.0
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.
15
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.
20 *
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.
27 *
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.
30 *
31 * For more info, check our website at http://www.gromacs.org
32 *
33 * And Hey:
34 * GROwing Monsters And Cloning Shrimps
35 */
36 #ifdef HAVE_CONFIG_H
37 #include <config.h>
38 #endif
39
40 #include <string.h>
41 #include "typedefs.h"
42 #include "string2.h"
43 #include "mdebin.h"
44 #include "smalloc.h"
45 #include "physics.h"
46 #include "enxio.h"
47 #include "vec.h"
48 #include "disre.h"
49 #include "main.h"
50 #include "network.h"
51 #include "names.h"
52 #include "orires.h"
53 #include "constr.h"
54 #include "mtop_util.h"
55 #include "xvgr.h"
56 #include "gmxfio.h"
57
58 static const char *conrmsd_nm[] = { "Constr. rmsd", "Constr.2 rmsd" };
59
60 static const char *boxs_nm[] = { "Box-X", "Box-Y", "Box-Z" };
61
62 static const char *tricl_boxs_nm[] = { "Box-XX", "Box-YX", "Box-YY",
63 "Box-ZX", "Box-ZY", "Box-ZZ" };
64
65 static const char *vol_nm[] = { "Volume" };
66
67 static const char *dens_nm[] = {"Density" };
68
69 static const char *pv_nm[] = {"pV" };
70
71 static const char *enthalpy_nm[] = {"Enthalpy" };
72
73 static const char *boxvel_nm[] = {
74 "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
75 "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
76 };
77
78 #define NBOXS asize(boxs_nm)
79 #define NTRICLBOXS asize(tricl_boxs_nm)
80
81 static bool bTricl,bDynBox;
82 static int f_nre=0,epc,etc,nCrmsd;
83
84 t_mdebin *init_mdebin(ener_file_t fp_ene,
85 const gmx_mtop_t *mtop,
86 const t_inputrec *ir)
87 {
88 const char *ener_nm[F_NRE];
89 static const char *vir_nm[] = {
90 "Vir-XX", "Vir-XY", "Vir-XZ",
91 "Vir-YX", "Vir-YY", "Vir-YZ",
92 "Vir-ZX", "Vir-ZY", "Vir-ZZ"
93 };
94 static const char *sv_nm[] = {
95 "ShakeVir-XX", "ShakeVir-XY", "ShakeVir-XZ",
96 "ShakeVir-YX", "ShakeVir-YY", "ShakeVir-YZ",
97 "ShakeVir-ZX", "ShakeVir-ZY", "ShakeVir-ZZ"
98 };
99 static const char *fv_nm[] = {
100 "ForceVir-XX", "ForceVir-XY", "ForceVir-XZ",
101 "ForceVir-YX", "ForceVir-YY", "ForceVir-YZ",
102 "ForceVir-ZX", "ForceVir-ZY", "ForceVir-ZZ"
103 };
104 static const char *pres_nm[] = {
105 "Pres-XX","Pres-XY","Pres-XZ",
106 "Pres-YX","Pres-YY","Pres-YZ",
107 "Pres-ZX","Pres-ZY","Pres-ZZ"
108 };
109 static const char *surft_nm[] = {
110 "#Surf*SurfTen"
111 };
112 static const char *mu_nm[] = {
113 "Mu-X", "Mu-Y", "Mu-Z"
114 };
115 static const char *vcos_nm[] = {
116 "2CosZ*Vel-X"
117 };
118 static const char *visc_nm[] = {
119 "1/Viscosity"
120 };
121 static const char *baro_nm[] = {
122 "Barostat"
123 };
124
125 char **grpnms;
126 const gmx_groups_t *groups;
127 char **gnm;
128 char buf[256];
129 const char *bufi;
130 t_mdebin *md;
131 int i,j,ni,nj,n,nh,k,kk,ncon,nset;
132 bool bBHAM,bNoseHoover,b14;
133
134 snew(md,1);
135
136 groups = &mtop->groups;
137
138 bBHAM = (mtop->ffparams.functype[0] == F_BHAM);
139 b14 = (gmx_mtop_ftype_count(mtop,F_LJ14) > 0 ||
140 gmx_mtop_ftype_count(mtop,F_LJC14_Q) > 0);
141
142 ncon = gmx_mtop_ftype_count(mtop,F_CONSTR);
143 nset = gmx_mtop_ftype_count(mtop,F_SETTLE);
144 md->bConstr = (ncon > 0 || nset > 0);
145 md->bConstrVir = FALSE;
146 if (md->bConstr) {
147 if (ncon > 0 && ir->eConstrAlg == econtLINCS) {
148 if (ir->eI == eiSD2)
149 md->nCrmsd = 2;
150 else
151 md->nCrmsd = 1;
152 }
153 md->bConstrVir = (getenv("GMX_CONSTRAINTVIR") != NULL);
154 } else {
155 md->nCrmsd = 0;
156 }
157
158 /* Energy monitoring */
159 for(i=0;i<egNR;i++)
160 {
161 md->bEInd[i]=FALSE;
162 }
163
164 for(i=0; i<F_NRE; i++) {
165 md->bEner[i] = FALSE;
166 if (i == F_LJ)
167 md->bEner[i] = !bBHAM;
168 else if (i == F_BHAM)
169 md->bEner[i] = bBHAM;
170 else if (i == F_EQM)
171 md->bEner[i] = ir->bQMMM;
172 else if (i == F_COUL_LR)
173 md->bEner[i] = (ir->rcoulomb > ir->rlist);
174 else if (i == F_LJ_LR)
175 md->bEner[i] = (!bBHAM && ir->rvdw > ir->rlist);
176 else if (i == F_BHAM_LR)
177 md->bEner[i] = (bBHAM && ir->rvdw > ir->rlist);
178 else if (i == F_RF_EXCL)
179 md->bEner[i] = (EEL_RF(ir->coulombtype) && ir->coulombtype != eelRF_NEC);
180 else if (i == F_COUL_RECIP)
181 md->bEner[i] = EEL_FULL(ir->coulombtype);
182 else if (i == F_LJ14)
183 md->bEner[i] = b14;
184 else if (i == F_COUL14)
185 md->bEner[i] = b14;
186 else if (i == F_LJC14_Q || i == F_LJC_PAIRS_NB)
187 md->bEner[i] = FALSE;
188 else if ((i == F_DVDL) || (i == F_DKDL))
189 md->bEner[i] = (ir->efep != efepNO);
190 else if (i == F_DHDL_CON)
191 md->bEner[i] = (ir->efep != efepNO && md->bConstr);
192 else if ((interaction_function[i].flags & IF_VSITE) ||
193 (i == F_CONSTR) || (i == F_CONSTRNC) || (i == F_SETTLE))
194 md->bEner[i] = FALSE;
195 else if ((i == F_COUL_SR) || (i == F_EPOT) || (i == F_PRES) || (i==F_EQM))
196 md->bEner[i] = TRUE;
197 else if ((i == F_ETOT) || (i == F_EKIN) || (i == F_TEMP))
198 md->bEner[i] = EI_DYNAMICS(ir->eI);
199 else if (i==F_VTEMP)
200 md->bEner[i] = (EI_DYNAMICS(ir->eI) && getenv("GMX_VIRIAL_TEMPERATURE"));
201 else if (i == F_DISPCORR || i == F_PDISPCORR)
202 md->bEner[i] = (ir->eDispCorr != edispcNO);
203 else if (i == F_DISRESVIOL)
204 md->bEner[i] = (gmx_mtop_ftype_count(mtop,F_DISRES) > 0);
205 else if (i == F_ORIRESDEV)
206 md->bEner[i] = (gmx_mtop_ftype_count(mtop,F_ORIRES) > 0);
207 else if (i == F_CONNBONDS)
208 md->bEner[i] = FALSE;
209 else if (i == F_COM_PULL)
210 md->bEner[i] = (ir->ePull == epullUMBRELLA || ir->ePull == epullCONST_F);
211 else if (i == F_ECONSERVED)
212 md->bEner[i] = ((ir->etc == etcNOSEHOOVER || ir->etc == etcVRESCALE) &&
213 (ir->epc == epcNO || ir->epc==epcMTTK));
214 else
215 md->bEner[i] = (gmx_mtop_ftype_count(mtop,i) > 0);
216 }
217
218 md->f_nre=0;
219 for(i=0; i<F_NRE; i++)
220 {
221 if (md->bEner[i])
222 {
223 /* FIXME: The constness should not be cast away */
224 /*ener_nm[f_nre]=(char *)interaction_function[i].longname;*/
225 ener_nm[md->f_nre]=interaction_function[i].longname;
226 md->f_nre++;
227 }
228 }
229
230 md->epc = ir->epc;
231 for (i=0;i<DIM;i++)
232 {
233 for (j=0;j<DIM;j++)
234 {
235 md->ref_p[i][j] = ir->ref_p[i][j];
236 }
237 }
238 md->bTricl = TRICLINIC(ir->compress) || TRICLINIC(ir->deform);
239 md->bDynBox = DYNAMIC_BOX(*ir);
240 md->etc = ir->etc;
241 md->bNHC_trotter = IR_NVT_TROTTER(ir);
242 md->bMTTK = IR_NPT_TROTTER(ir);
243
244 md->ebin = mk_ebin();
245 /* Pass NULL for unit to let get_ebin_space determine the units
246 * for interaction_function[i].longname
247 */
248 md->ie = get_ebin_space(md->ebin,md->f_nre,ener_nm,NULL);
249 if (md->nCrmsd)
250 {
251 /* This should be called directly after the call for md->ie,
252 * such that md->iconrmsd follows directly in the list.
253 */
254 md->iconrmsd = get_ebin_space(md->ebin,md->nCrmsd,conrmsd_nm,"");
255 }
256 if (md->bDynBox)
257 {
258 md->ib = get_ebin_space(md->ebin, md->bTricl ? NTRICLBOXS :
259 NBOXS, md->bTricl ? tricl_boxs_nm : boxs_nm,
260 unit_length);
261 md->ivol = get_ebin_space(md->ebin, 1, vol_nm, unit_volume);
262 md->idens = get_ebin_space(md->ebin, 1, dens_nm, unit_density_SI);
263 md->ipv = get_ebin_space(md->ebin, 1, pv_nm, unit_energy);
264 md->ienthalpy = get_ebin_space(md->ebin, 1, enthalpy_nm, unit_energy);
265 }
266 if (md->bConstrVir)
267 {
268 md->isvir = get_ebin_space(md->ebin,asize(sv_nm),sv_nm,unit_energy);
269 md->ifvir = get_ebin_space(md->ebin,asize(fv_nm),fv_nm,unit_energy);
270 }
271 md->ivir = get_ebin_space(md->ebin,asize(vir_nm),vir_nm,unit_energy);
272 md->ipres = get_ebin_space(md->ebin,asize(pres_nm),pres_nm,unit_pres_bar);
273 md->isurft = get_ebin_space(md->ebin,asize(surft_nm),surft_nm,
274 unit_surft_bar);
275 if (md->epc == epcPARRINELLORAHMAN || md->epc == epcMTTK)
276 {
277 md->ipc = get_ebin_space(md->ebin,md->bTricl ? 6 : 3,
278 boxvel_nm,unit_vel);
279 }
280 md->imu = get_ebin_space(md->ebin,asize(mu_nm),mu_nm,unit_dipole_D);
281 if (ir->cos_accel != 0)
282 {
283 md->ivcos = get_ebin_space(md->ebin,asize(vcos_nm),vcos_nm,unit_vel);
284 md->ivisc = get_ebin_space(md->ebin,asize(visc_nm),visc_nm,
285 unit_invvisc_SI);
286 }
287
288 /* Energy monitoring */
289 for(i=0;i<egNR;i++)
290 {
291 md->bEInd[i] = FALSE;
292 }
293 md->bEInd[egCOULSR] = TRUE;
294 md->bEInd[egLJSR ] = TRUE;
295
296 if (ir->rcoulomb > ir->rlist)
297 {
298 md->bEInd[egCOULLR] = TRUE;
299 }
300 if (!bBHAM)
301 {
302 if (ir->rvdw > ir->rlist)
303 {
304 md->bEInd[egLJLR] = TRUE;
305 }
306 }
307 else
308 {
309 md->bEInd[egLJSR] = FALSE;
310 md->bEInd[egBHAMSR] = TRUE;
311 if (ir->rvdw > ir->rlist)
312 {
313 md->bEInd[egBHAMLR] = TRUE;
314 }
315 }
316 if (b14)
317 {
318 md->bEInd[egLJ14] = TRUE;
319 md->bEInd[egCOUL14] = TRUE;
320 }
321 md->nEc=0;
322 for(i=0; (i<egNR); i++)
323 {
324 if (md->bEInd[i])
325 {
326 md->nEc++;
327 }
328 }
329
330 n=groups->grps[egcENER].nr;
331 md->nEg=n;
332 md->nE=(n*(n+1))/2;
333 snew(md->igrp,md->nE);
334 if (md->nE > 1)
335 {
336 n=0;
337 snew(gnm,md->nEc);
338 for(k=0; (k<md->nEc); k++)
339 {
340 snew(gnm[k],STRLEN);
341 }
342 for(i=0; (i<groups->grps[egcENER].nr); i++)
343 {
344 ni=groups->grps[egcENER].nm_ind[i];
345 for(j=i; (j<groups->grps[egcENER].nr); j++)
346 {
347 nj=groups->grps[egcENER].nm_ind[j];
348 for(k=kk=0; (k<egNR); k++)
349 {
350 if (md->bEInd[k])
351 {
352 sprintf(gnm[kk],"%s:%s-%s",egrp_nm[k],
353 *(groups->grpname[ni]),*(groups->grpname[nj]));
354 kk++;
355 }
356 }
357 md->igrp[n]=get_ebin_space(md->ebin,md->nEc,
358 (const char **)gnm,unit_energy);
359 n++;
360 }
361 }
362 for(k=0; (k<md->nEc); k++)
363 {
364 sfree(gnm[k]);
365 }
366 sfree(gnm);
367
368 if (n != md->nE)
369 {
370 gmx_incons("Number of energy terms wrong");
371 }
372 }
373
374
375 md->nTC=groups->grps[egcTC].nr;
376 md->nNHC = ir->opts.nhchainlength; /* shorthand for number of NH chains */
377 if (md->bMTTK)
378 {
379 md->nTCP = 1; /* assume only one possible coupling system for barostat for now */
380 }
381 else
382 {
383 md->nTCP = 0;
384 }
385
386 if (md->etc == etcNOSEHOOVER) {
387 if (md->bNHC_trotter) {
388 md->mde_n = 2*md->nNHC*md->nTC;
389 }
390 else
391 {
392 md->mde_n = 2*md->nTC;
393 }
394 if (md->epc == epcMTTK)
395 {
396 md->mdeb_n = 2*md->nNHC*md->nTCP;
397 }
398 } else {
399 md->mde_n = md->nTC;
400 md->mdeb_n = 0;
401 }
402
403 snew(md->tmp_r,md->mde_n);
404 snew(md->tmp_v,md->mde_n);
405 snew(md->grpnms,md->mde_n);
406 grpnms = md->grpnms;
407
408 for(i=0; (i<md->nTC); i++)
409 {
410 ni=groups->grps[egcTC].nm_ind[i];
411 sprintf(buf,"T-%s",*(groups->grpname[ni]));
412 grpnms[i]=strdup(buf);
413 }
414 md->itemp=get_ebin_space(md->ebin,md->nTC,(const char **)grpnms,
415 unit_temp_K);
416
417 bNoseHoover = (getenv("GMX_NOSEHOOVER_CHAINS") != NULL); /* whether to print Nose-Hoover chains */
418
419 if (md->etc == etcNOSEHOOVER)
420 {
421 if (bNoseHoover)
422 {
423 if (md->bNHC_trotter)
424 {
425 for(i=0; (i<md->nTC); i++)
426 {
427 ni=groups->grps[egcTC].nm_ind[i];
428 bufi = *(groups->grpname[ni]);
429 for(j=0; (j<md->nNHC); j++)
430 {
431 sprintf(buf,"Xi-%d-%s",j,bufi);
432 grpnms[2*(i*md->nNHC+j)]=strdup(buf);
433 sprintf(buf,"vXi-%d-%s",j,bufi);
434 grpnms[2*(i*md->nNHC+j)+1]=strdup(buf);
435 }
436 }
437 md->itc=get_ebin_space(md->ebin,md->mde_n,(const char **)grpnms,unit_invtime);
438 if (md->bMTTK)
439 {
440 for(i=0; (i<md->nTCP); i++)
441 {
442 bufi = baro_nm[0]; /* All barostat DOF's together for now. */
443 for(j=0; (j<md->nNHC); j++)
444 {
445 sprintf(buf,"Xi-%d-%s",j,bufi);
446 grpnms[2*(i*md->nNHC+j)]=strdup(buf);
447 sprintf(buf,"vXi-%d-%s",j,bufi);
448 grpnms[2*(i*md->nNHC+j)+1]=strdup(buf);
449 }
450 }
451 md->itcb=get_ebin_space(md->ebin,md->mdeb_n,(const char **)grpnms,unit_invtime);
452 }
453 }
454 else
455 {
456 for(i=0; (i<md->nTC); i++)
457 {
458 ni=groups->grps[egcTC].nm_ind[i];
459 bufi = *(groups->grpname[ni]);
460 sprintf(buf,"Xi-%s",bufi);
461 grpnms[2*i]=strdup(buf);
462 sprintf(buf,"vXi-%s",bufi);
463 grpnms[2*i+1]=strdup(buf);
464 }
465 md->itc=get_ebin_space(md->ebin,md->mde_n,(const char **)grpnms,unit_invtime);
466 }
467 }
468 }
469 else if (md->etc == etcBERENDSEN || md->etc == etcYES ||
470 md->etc == etcVRESCALE)
471 {
472 for(i=0; (i<md->nTC); i++)
473 {
474 ni=groups->grps[egcTC].nm_ind[i];
475 sprintf(buf,"Lamb-%s",*(groups->grpname[ni]));
476 grpnms[i]=strdup(buf);
477 }
478 md->itc=get_ebin_space(md->ebin,md->mde_n,(const char **)grpnms,"");
479 }
480
481 sfree(grpnms);
482
483
484 md->nU=groups->grps[egcACC].nr;
485 if (md->nU > 1)
486 {
487 snew(grpnms,3*md->nU);
488 for(i=0; (i<md->nU); i++)
489 {
490 ni=groups->grps[egcACC].nm_ind[i];
491 sprintf(buf,"Ux-%s",*(groups->grpname[ni]));
492 grpnms[3*i+XX]=strdup(buf);
493 sprintf(buf,"Uy-%s",*(groups->grpname[ni]));
494 grpnms[3*i+YY]=strdup(buf);
495 sprintf(buf,"Uz-%s",*(groups->grpname[ni]));
496 grpnms[3*i+ZZ]=strdup(buf);
497 }
498 md->iu=get_ebin_space(md->ebin,3*md->nU,(const char **)grpnms,unit_vel);
499 sfree(grpnms);
500 }
501
502 if ( fp_ene )
503 {
504 do_enxnms(fp_ene,&md->ebin->nener,&md->ebin->enm);
505 }
506
507 md->print_grpnms=NULL;
508
509 return md;
510 }
511
512 FILE *open_dhdl(const char *filename,const t_inputrec *ir,
513 const output_env_t oenv)
514 {
515 FILE *fp;
516 const char *dhdl="dH/d\\lambda",*deltag="\\DeltaH",*lambda="\\lambda";
517 char title[STRLEN],label_x[STRLEN],label_y[STRLEN];
518 int nsets,s;
519 char **setname,buf[STRLEN];
520
521 sprintf(label_x,"%s (%s)","Time",unit_time);
522 if (ir->n_flambda == 0)
523 {
524 sprintf(title,"%s",dhdl);
525 sprintf(label_y,"%s (%s %s)",
526 dhdl,unit_energy,"[\\lambda]\\S-1\\N");
527 }
528 else
529 {
530 sprintf(title,"%s, %s",dhdl,deltag);
531 sprintf(label_y,"(%s)",unit_energy);
532 }
533 fp = xvgropen(filename,title,label_x,label_y,oenv);
534
535 sprintf(buf,"T = %g (K)",ir->opts.ref_t[0]);
536 xvgr_subtitle(fp,buf,oenv);
537
538 if (ir->n_flambda > 0)
539 {
540 /* g_bar has to determine the lambda values used in this simulation
541 * from this xvg legend.
542 */
543 nsets = 1 + ir->n_flambda;
544 snew(setname,nsets);
545 sprintf(buf,"%s %s %g",dhdl,lambda,ir->init_lambda);
546 setname[0] = strdup(buf);
547 for(s=1; s<nsets; s++)
548 {
549 sprintf(buf,"%s %s %g",deltag,lambda,ir->flambda[s-1]);
550 setname[s] = strdup(buf);
551 }
552 xvgr_legend(fp,nsets,setname,oenv);
553
554 for(s=0; s<nsets; s++)
555 {
556 sfree(setname[s]);
557 }
558 sfree(setname);
559 }
560
561 return fp;
562 }
563
564 static void copy_energy(t_mdebin *md, real e[],real ecpy[])
565 {
566 int i,j;
567
568 for(i=j=0; (i<F_NRE); i++)
569 if (md->bEner[i])
570 ecpy[j++] = e[i];
571 if (j != md->f_nre)
572 gmx_incons("Number of energy terms wrong");
573 }
574
575 void upd_mdebin(t_mdebin *md,FILE *fp_dhdl,
576 bool bSum,
577 double time,
578 real tmass,
579 gmx_enerdata_t *enerd,
580 t_state *state,
581 matrix box,
582 tensor svir,
583 tensor fvir,
584 tensor vir,
585 tensor pres,
586 gmx_ekindata_t *ekind,
587 rvec mu_tot,
588 gmx_constr_t constr)
589 {
590 int i,j,k,kk,m,n,gid;
591 real crmsd[2],tmp6[6];
592 real bs[NTRICLBOXS],vol,dens,pv,enthalpy;
593 real eee[egNR];
594 real ecopy[F_NRE];
595 real tmp;
596 bool bNoseHoover;
597
598 /* Do NOT use the box in the state variable, but the separate box provided
599 * as an argument. This is because we sometimes need to write the box from
600 * the last timestep to match the trajectory frames.
601 */
602 copy_energy(md, enerd->term,ecopy);
603 add_ebin(md->ebin,md->ie,md->f_nre,ecopy,bSum);
604 if (md->nCrmsd)
605 {
606 crmsd[0] = constr_rmsd(constr,FALSE);
607 if (md->nCrmsd > 1)
608 {
609 crmsd[1] = constr_rmsd(constr,TRUE);
610 }
611 add_ebin(md->ebin,md->iconrmsd,md->nCrmsd,crmsd,FALSE);
612 }
613 if (md->bDynBox)
614 {
615 if(md->bTricl)
616 {
617 bs[0] = box[XX][XX];
618 bs[1] = box[YY][XX];
619 bs[2] = box[YY][YY];
620 bs[3] = box[ZZ][XX];
621 bs[4] = box[ZZ][YY];
622 bs[5] = box[ZZ][ZZ];
623 }
624 else
625 {
626 bs[0] = box[XX][XX];
627 bs[1] = box[YY][YY];
628 bs[2] = box[ZZ][ZZ];
629 }
630 vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
631 dens = (tmass*AMU)/(vol*NANO*NANO*NANO);
632
633 /* This is pV (in kJ/mol). The pressure is the reference pressure,
634 not the instantaneous pressure */
635 pv = 0;
636 for (i=0;i<DIM;i++)
637 {
638 for (j=0;j<DIM;j++)
639 {
640 if (i>j)
641 {
642 pv += box[i][j]*md->ref_p[i][j]/PRESFAC;
643 }
644 else
645 {
646 pv += box[j][i]*md->ref_p[j][i]/PRESFAC;
647 }
648 }
649 }
650 add_ebin(md->ebin,md->ib ,NBOXS,bs ,bSum);
651 add_ebin(md->ebin,md->ivol ,1 ,&vol ,bSum);
652 add_ebin(md->ebin,md->idens,1 ,&dens,bSum);
653 add_ebin(md->ebin,md->ipv ,1 ,&pv ,bSum);
654 enthalpy = pv + enerd->term[F_ETOT];
655 add_ebin(md->ebin,md->ienthalpy ,1 ,&enthalpy ,bSum);
656 }
657 if (md->bConstrVir)
658 {
659 add_ebin(md->ebin,md->isvir,9,svir[0],bSum);
660 add_ebin(md->ebin,md->ifvir,9,fvir[0],bSum);
661 }
662 add_ebin(md->ebin,md->ivir,9,vir[0],bSum);
663 add_ebin(md->ebin,md->ipres,9,pres[0],bSum);
664 tmp = (pres[ZZ][ZZ]-(pres[XX][XX]+pres[YY][YY])*0.5)*box[ZZ][ZZ];
665 add_ebin(md->ebin,md->isurft,1,&tmp,bSum);
666 if (md->epc == epcPARRINELLORAHMAN || md->epc == epcMTTK)
667 {
668 tmp6[0] = state->boxv[XX][XX];
669 tmp6[1] = state->boxv[YY][YY];
670 tmp6[2] = state->boxv[ZZ][ZZ];
671 tmp6[3] = state->boxv[YY][XX];
672 tmp6[4] = state->boxv[ZZ][XX];
673 tmp6[5] = state->boxv[ZZ][YY];
674 add_ebin(md->ebin,md->ipc,md->bTricl ? 6 : 3,tmp6,bSum);
675 }
676 add_ebin(md->ebin,md->imu,3,mu_tot,bSum);
677 if (ekind && ekind->cosacc.cos_accel != 0)
678 {
679 vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
680 dens = (tmass*AMU)/(vol*NANO*NANO*NANO);
681 add_ebin(md->ebin,md->ivcos,1,&(ekind->cosacc.vcos),bSum);
682 /* 1/viscosity, unit 1/(kg m^-1 s^-1) */
683 tmp = 1/(ekind->cosacc.cos_accel/(ekind->cosacc.vcos*PICO)
684 *vol*sqr(box[ZZ][ZZ]*NANO/(2*M_PI)));
685 add_ebin(md->ebin,md->ivisc,1,&tmp,bSum);
686 }
687 if (md->nE > 1)
688 {
689 n=0;
690 for(i=0; (i<md->nEg); i++)
691 {
692 for(j=i; (j<md->nEg); j++)
693 {
694 gid=GID(i,j,md->nEg);
695 for(k=kk=0; (k<egNR); k++)
696 {
697 if (md->bEInd[k])
698 {
699 eee[kk++] = enerd->grpp.ener[k][gid];
700 }
701 }
702 add_ebin(md->ebin,md->igrp[n],md->nEc,eee,bSum);
703 n++;
704 }
705 }
706 }
707
708 if (ekind)
709 {
710 for(i=0; (i<md->nTC); i++)
711 {
712 md->tmp_r[i] = ekind->tcstat[i].T;
713 }
714 add_ebin(md->ebin,md->itemp,md->nTC,md->tmp_r,bSum);
715
716 bNoseHoover = (getenv("GMX_NOSEHOOVER_CHAINS") != NULL); /* whether to print Nose-Hoover chains */
717
718 if (md->etc == etcNOSEHOOVER)
719 {
720 if (bNoseHoover)
721 {
722 if (md->bNHC_trotter)
723 {
724 for(i=0; (i<md->nTC); i++)
725 {
726 for (j=0;j<md->nNHC;j++)
727 {
728 k = i*md->nNHC+j;
729 md->tmp_r[2*k] = state->nosehoover_xi[k];
730 md->tmp_r[2*k+1] = state->nosehoover_vxi[k];
731 }
732 }
733 add_ebin(md->ebin,md->itc,md->mde_n,md->tmp_r,bSum);
734
735 if (md->bMTTK) {
736 for(i=0; (i<md->nTCP); i++)
737 {
738 for (j=0;j<md->nNHC;j++)
739 {
740 k = i*md->nNHC+j;
741 md->tmp_r[2*k] = state->nhpres_xi[k];
742 md->tmp_r[2*k+1] = state->nhpres_vxi[k];
743 }
744 }
745 add_ebin(md->ebin,md->itcb,md->mdeb_n,md->tmp_r,bSum);
746 }
747
748 }
749 else
750 {
751 for(i=0; (i<md->nTC); i++)
752 {
753 md->tmp_r[2*i] = state->nosehoover_xi[i];
754 md->tmp_r[2*i+1] = state->nosehoover_vxi[i];
755 }
756 add_ebin(md->ebin,md->itc,md->mde_n,md->tmp_r,bSum);
757 }
758 }
759 }
760 else if (md->etc == etcBERENDSEN || md->etc == etcYES || md->etc == etcVRESCALE)
761 {
762 for(i=0; (i<md->nTC); i++)
763 {
764 md->tmp_r[i] = ekind->tcstat[i].lambda;
765 }
766 add_ebin(md->ebin,md->itc,md->nTC,md->tmp_r,bSum);
767 }
768 }
769
770 if (ekind && md->nU > 1)
771 {
772 for(i=0; (i<md->nU); i++)
773 {
774 copy_rvec(ekind->grpstat[i].u,md->tmp_v[i]);
775 }
776 add_ebin(md->ebin,md->iu,3*md->nU,md->tmp_v[0],bSum);
777 }
778
779 ebin_increase_count(md->ebin,bSum);
780
781 if (fp_dhdl)
782 {
783 fprintf(fp_dhdl,"%.4f %g",
784 time,
785 enerd->term[F_DVDL]+enerd->term[F_DKDL]+enerd->term[F_DHDL_CON]);
786 for(i=1; i<enerd->n_lambda; i++)
787 {
788 fprintf(fp_dhdl," %g",
789 enerd->enerpart_lambda[i]-enerd->enerpart_lambda[0]);
790 }
791 fprintf(fp_dhdl,"\n");
792 }
793 }
794
795 void upd_mdebin_step(t_mdebin *md)
796 {
797 ebin_increase_count(md->ebin,FALSE);
798 }
799
800 static void npr(FILE *log,int n,char c)
801 {
802 for(; (n>0); n--) fprintf(log,"%c",c);
803 }
804
805 static void pprint(FILE *log,const char *s,t_mdebin *md)
806 {
807 char CHAR='#';
808 int slen;
809 char buf1[22],buf2[22];
810
811 slen = strlen(s);
812 fprintf(log,"\t<====== ");
813 npr(log,slen,CHAR);
814 fprintf(log," ==>\n");
815 fprintf(log,"\t<==== %s ====>\n",s);
816 fprintf(log,"\t<== ");
817 npr(log,slen,CHAR);
818 fprintf(log," ======>\n\n");
819
820 fprintf(log,"\tStatistics over %s steps using %s frames\n",
821 gmx_step_str(md->ebin->nsteps_sim,buf1),
822 gmx_step_str(md->ebin->nsum_sim,buf2));
823 fprintf(log,"\n");
824 }
825
826 void print_ebin_header(FILE *log,gmx_large_int_t steps,double time,real lamb)
827 {
828 char buf[22];
829
830 fprintf(log," %12s %12s %12s\n"
831 " %12s %12.5f %12.5f\n\n",
832 "Step","Time","Lambda",gmx_step_str(steps,buf),time,lamb);
833 }
834
835 void print_ebin(ener_file_t fp_ene,bool bEne,bool bDR,bool bOR,
836 FILE *log,
837 gmx_large_int_t step,double time,
838 int mode,bool bCompact,
839 t_mdebin *md,t_fcdata *fcd,
840 gmx_groups_t *groups,t_grpopts *opts)
841 {
842 /*static char **grpnms=NULL;*/
843 char buf[246];
844 int i,j,n,ni,nj,ndr,nor;
845 int nr[enxNR];
846 real *block[enxNR];
847 t_enxframe fr;
848
849 switch (mode)
850 {
851 case eprNORMAL:
852 fr.t = time;
853 fr.step = step;
854 fr.nsteps = md->ebin->nsteps;
855 fr.nsum = md->ebin->nsum;
856 fr.nre = (bEne) ? md->ebin->nener : 0;
857 fr.ener = md->ebin->e;
858 fr.ndisre = bDR ? fcd->disres.npair : 0;
859 fr.disre_rm3tav = fcd->disres.rm3tav;
860 fr.disre_rt = fcd->disres.rt;
861 /* Optional additional blocks */
862 for(i=0; i<enxNR; i++)
863 {
864 nr[i] = 0;
865 }
866 if (fcd->orires.nr > 0 && bOR)
867 {
868 diagonalize_orires_tensors(&(fcd->orires));
869 nr[enxOR] = fcd->orires.nr;
870 block[enxOR] = fcd->orires.otav;
871 nr[enxORI] = (fcd->orires.oinsl != fcd->orires.otav) ?
872 fcd->orires.nr : 0;
873 block[enxORI] = fcd->orires.oinsl;
874 nr[enxORT] = fcd->orires.nex*12;
875 block[enxORT] = fcd->orires.eig;
876 }
877 fr.nblock = 0;
878 for(i=0; i<enxNR; i++)
879 {
880 if (nr[i] > 0)
881 {
882 fr.nblock = i + 1;
883 }
884 }
885 fr.nr = nr;
886 fr.block = block;
887 if (fr.nre || fr.ndisre || fr.nr[enxOR] || fr.nr[enxORI])
888 {
889 do_enx(fp_ene,&fr);
890 gmx_fio_check_file_position(enx_file_pointer(fp_ene));
891 if (fr.nre)
892 {
893 /* We have stored the sums, so reset the sum history */
894 reset_ebin_sums(md->ebin);
895 }
896 }
897 break;
898 case eprAVER:
899 if (log)
900 {
901 pprint(log,"A V E R A G E S",md);
902 }
903 break;
904 case eprRMS:
905 if (log)
906 {
907 pprint(log,"R M S - F L U C T U A T I O N S",md);
908 }
909 break;
910 default:
911 gmx_fatal(FARGS,"Invalid print mode (%d)",mode);
912 }
913
914 if (log)
915 {
916 for(i=0;i<opts->ngtc;i++)
917 {
918 if(opts->annealing[i]!=eannNO)
919 {
920 fprintf(log,"Current ref_t for group %s: %8.1f\n",
921 *(groups->grpname[groups->grps[egcTC].nm_ind[i]]),
922 opts->ref_t[i]);
923 }
924 }
925 if (mode==eprNORMAL && fcd->orires.nr>0)
926 {
927 print_orires_log(log,&(fcd->orires));
928 }
929 fprintf(log," Energies (%s)\n",unit_energy);
930 pr_ebin(log,md->ebin,md->ie,md->f_nre+md->nCrmsd,5,mode,TRUE);
931 fprintf(log,"\n");
932
933 if (!bCompact)
934 {
935 if (md->bDynBox)
936 {
937 pr_ebin(log,md->ebin,md->ib, md->bTricl ? NTRICLBOXS : NBOXS,5,
938 mode,TRUE);
939 fprintf(log,"\n");
940 }
941 if (md->bConstrVir)
942 {
943 fprintf(log," Constraint Virial (%s)\n",unit_energy);
944 pr_ebin(log,md->ebin,md->isvir,9,3,mode,FALSE);
945 fprintf(log,"\n");
946 fprintf(log," Force Virial (%s)\n",unit_energy);
947 pr_ebin(log,md->ebin,md->ifvir,9,3,mode,FALSE);
948 fprintf(log,"\n");
949 }
950 fprintf(log," Total Virial (%s)\n",unit_energy);
951 pr_ebin(log,md->ebin,md->ivir,9,3,mode,FALSE);
952 fprintf(log,"\n");
953 fprintf(log," Pressure (%s)\n",unit_pres_bar);
954 pr_ebin(log,md->ebin,md->ipres,9,3,mode,FALSE);
955 fprintf(log,"\n");
956 fprintf(log," Total Dipole (%s)\n",unit_dipole_D);
957 pr_ebin(log,md->ebin,md->imu,3,3,mode,FALSE);
958 fprintf(log,"\n");
959
960 if (md->nE > 1)
961 {
962 if (md->print_grpnms==NULL)
963 {
964 snew(md->print_grpnms,md->nE);
965 n=0;
966 for(i=0; (i<md->nEg); i++)
967 {
968 ni=groups->grps[egcENER].nm_ind[i];
969 for(j=i; (j<md->nEg); j++)
970 {
971 nj=groups->grps[egcENER].nm_ind[j];
972 sprintf(buf,"%s-%s",*(groups->grpname[ni]),
973 *(groups->grpname[nj]));
974 md->print_grpnms[n++]=strdup(buf);
975 }
976 }
977 }
978 sprintf(buf,"Epot (%s)",unit_energy);
979 fprintf(log,"%15s ",buf);
980 for(i=0; (i<egNR); i++)
981 {
982 if (md->bEInd[i])
983 {
984 fprintf(log,"%12s ",egrp_nm[i]);
985 }
986 }
987 fprintf(log,"\n");
988 for(i=0; (i<md->nE); i++)
989 {
990 fprintf(log,"%15s",md->print_grpnms[i]);
991 pr_ebin(log,md->ebin,md->igrp[i],md->nEc,md->nEc,mode,
992 FALSE);
993 }
994 fprintf(log,"\n");
995 }
996 if (md->nTC > 1)
997 {
998 pr_ebin(log,md->ebin,md->itemp,md->nTC,4,mode,TRUE);
999 fprintf(log,"\n");
1000 }
1001 if (md->nU > 1)
1002 {
1003 fprintf(log,"%15s %12s %12s %12s\n",
1004 "Group","Ux","Uy","Uz");
1005 for(i=0; (i<md->nU); i++)
1006 {
1007 ni=groups->grps[egcACC].nm_ind[i];
1008 fprintf(log,"%15s",*groups->grpname[ni]);
1009 pr_ebin(log,md->ebin,md->iu+3*i,3,3,mode,FALSE);
1010 }
1011 fprintf(log,"\n");
1012 }
1013 }
1014
1015 }
1016 }
1017
1018 void
1019 init_energyhistory(energyhistory_t * enerhist)
1020 {
1021 enerhist->nener = 0;
1022
1023 enerhist->ener_ave = NULL;
1024 enerhist->ener_sum = NULL;
1025 enerhist->ener_sum_sim = NULL;
1026
1027 enerhist->nsteps = 0;
1028 enerhist->nsum = 0;
1029 enerhist->nsteps_sim = 0;
1030 enerhist->nsum_sim = 0;
1031 }
1032
1033 void
1034 update_energyhistory(energyhistory_t * enerhist,t_mdebin * mdebin)
1035 {
1036 int i;
1037
1038 enerhist->nsteps = mdebin->ebin->nsteps;
1039 enerhist->nsum = mdebin->ebin->nsum;
1040 enerhist->nsteps_sim = mdebin->ebin->nsteps_sim;
1041 enerhist->nsum_sim = mdebin->ebin->nsum_sim;
1042 enerhist->nener = mdebin->ebin->nener;
1043
1044 if (mdebin->ebin->nsum > 0)
1045 {
1046 /* Check if we need to allocate first */
1047 if(enerhist->ener_ave == NULL)
1048 {
1049 snew(enerhist->ener_ave,enerhist->nener);
1050 snew(enerhist->ener_sum,enerhist->nener);
1051 }
1052
1053 for(i=0;i<enerhist->nener;i++)
1054 {
1055 enerhist->ener_ave[i] = mdebin->ebin->e[i].eav;
1056 enerhist->ener_sum[i] = mdebin->ebin->e[i].esum;
1057 }
1058 }
1059
1060 if (mdebin->ebin->nsum_sim > 0)
1061 {
1062 /* Check if we need to allocate first */
1063 if(enerhist->ener_sum_sim == NULL)
1064 {
1065 snew(enerhist->ener_sum_sim,enerhist->nener);
1066 }
1067
1068 for(i=0;i<enerhist->nener;i++)
1069 {
1070 enerhist->ener_sum_sim[i] = mdebin->ebin->e_sim[i].esum;
1071 }
1072 }
1073 }
1074
1075 void
1076 restore_energyhistory_from_state(t_mdebin * mdebin,energyhistory_t * enerhist)
1077 {
1078 int i;
1079
1080 if ((enerhist->nsum > 0 || enerhist->nsum_sim > 0) &&
1081 mdebin->ebin->nener != enerhist->nener)
1082 {
1083 gmx_fatal(FARGS,"Mismatch between number of energies in run input (%d) and checkpoint file (%d).",
1084 mdebin->ebin->nener,enerhist->nener);
1085 }
1086
1087 mdebin->ebin->nsteps = enerhist->nsteps;
1088 mdebin->ebin->nsum = enerhist->nsum;
1089 mdebin->ebin->nsteps_sim = enerhist->nsteps_sim;
1090 mdebin->ebin->nsum_sim = enerhist->nsum_sim;
1091
1092 for(i=0; i<mdebin->ebin->nener; i++)
1093 {
1094 mdebin->ebin->e[i].eav =
1095 (enerhist->nsum > 0 ? enerhist->ener_ave[i] : 0);
1096 mdebin->ebin->e[i].esum =
1097 (enerhist->nsum > 0 ? enerhist->ener_sum[i] : 0);
1098 mdebin->ebin->e_sim[i].esum =
1099 (enerhist->nsum_sim > 0 ? enerhist->ener_sum_sim[i] : 0);
1100 }
1101 }
Atempt to write QMMM interface for Gamess-US