?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
minor fixes in ditribution files
[gromacs/qmmm-gamess-us.git] / src / mdlib / mdebin.c
blob528b6e9aa63f87bf6f6810b3ecfa0113ccb6ebe1
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 = gmx_fio_fopen(filename,"w+");
534 xvgr_header(fp,title,label_x,label_y,exvggtXNY,oenv);
535
536 if (ir->delta_lambda == 0)
537 {
538 sprintf(buf,"T = %g (K), %s = %g",
539 ir->opts.ref_t[0],lambda,ir->init_lambda);
540 }
541 else
542 {
543 sprintf(buf,"T = %g (K)",
544 ir->opts.ref_t[0]);
545 }
546 xvgr_subtitle(fp,buf,oenv);
547
548 if (ir->n_flambda > 0)
549 {
550 /* g_bar has to determine the lambda values used in this simulation
551 * from this xvg legend.
552 */
553 nsets = 1 + ir->n_flambda;
554 snew(setname,nsets);
555 sprintf(buf,"%s %s %g",dhdl,lambda,ir->init_lambda);
556 setname[0] = strdup(buf);
557 for(s=1; s<nsets; s++)
558 {
559 sprintf(buf,"%s %s %g",deltag,lambda,ir->flambda[s-1]);
560 setname[s] = strdup(buf);
561 }
562 xvgr_legend(fp,nsets,setname,oenv);
563
564 for(s=0; s<nsets; s++)
565 {
566 sfree(setname[s]);
567 }
568 sfree(setname);
569 }
570
571 return fp;
572 }
573
574 static void copy_energy(t_mdebin *md, real e[],real ecpy[])
575 {
576 int i,j;
577
578 for(i=j=0; (i<F_NRE); i++)
579 if (md->bEner[i])
580 ecpy[j++] = e[i];
581 if (j != md->f_nre)
582 gmx_incons("Number of energy terms wrong");
583 }
584
585 void upd_mdebin(t_mdebin *md,FILE *fp_dhdl,
586 bool bSum,
587 double time,
588 real tmass,
589 gmx_enerdata_t *enerd,
590 t_state *state,
591 matrix box,
592 tensor svir,
593 tensor fvir,
594 tensor vir,
595 tensor pres,
596 gmx_ekindata_t *ekind,
597 rvec mu_tot,
598 gmx_constr_t constr)
599 {
600 int i,j,k,kk,m,n,gid;
601 real crmsd[2],tmp6[6];
602 real bs[NTRICLBOXS],vol,dens,pv,enthalpy;
603 real eee[egNR];
604 real ecopy[F_NRE];
605 real tmp;
606 bool bNoseHoover;
607
608 /* Do NOT use the box in the state variable, but the separate box provided
609 * as an argument. This is because we sometimes need to write the box from
610 * the last timestep to match the trajectory frames.
611 */
612 copy_energy(md, enerd->term,ecopy);
613 add_ebin(md->ebin,md->ie,md->f_nre,ecopy,bSum);
614 if (md->nCrmsd)
615 {
616 crmsd[0] = constr_rmsd(constr,FALSE);
617 if (md->nCrmsd > 1)
618 {
619 crmsd[1] = constr_rmsd(constr,TRUE);
620 }
621 add_ebin(md->ebin,md->iconrmsd,md->nCrmsd,crmsd,FALSE);
622 }
623 if (md->bDynBox)
624 {
625 if(md->bTricl)
626 {
627 bs[0] = box[XX][XX];
628 bs[1] = box[YY][XX];
629 bs[2] = box[YY][YY];
630 bs[3] = box[ZZ][XX];
631 bs[4] = box[ZZ][YY];
632 bs[5] = box[ZZ][ZZ];
633 }
634 else
635 {
636 bs[0] = box[XX][XX];
637 bs[1] = box[YY][YY];
638 bs[2] = box[ZZ][ZZ];
639 }
640 vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
641 dens = (tmass*AMU)/(vol*NANO*NANO*NANO);
642
643 /* This is pV (in kJ/mol). The pressure is the reference pressure,
644 not the instantaneous pressure */
645 pv = 0;
646 for (i=0;i<DIM;i++)
647 {
648 for (j=0;j<DIM;j++)
649 {
650 if (i>j)
651 {
652 pv += box[i][j]*md->ref_p[i][j]/PRESFAC;
653 }
654 else
655 {
656 pv += box[j][i]*md->ref_p[j][i]/PRESFAC;
657 }
658 }
659 }
660 add_ebin(md->ebin,md->ib ,NBOXS,bs ,bSum);
661 add_ebin(md->ebin,md->ivol ,1 ,&vol ,bSum);
662 add_ebin(md->ebin,md->idens,1 ,&dens,bSum);
663 add_ebin(md->ebin,md->ipv ,1 ,&pv ,bSum);
664 enthalpy = pv + enerd->term[F_ETOT];
665 add_ebin(md->ebin,md->ienthalpy ,1 ,&enthalpy ,bSum);
666 }
667 if (md->bConstrVir)
668 {
669 add_ebin(md->ebin,md->isvir,9,svir[0],bSum);
670 add_ebin(md->ebin,md->ifvir,9,fvir[0],bSum);
671 }
672 add_ebin(md->ebin,md->ivir,9,vir[0],bSum);
673 add_ebin(md->ebin,md->ipres,9,pres[0],bSum);
674 tmp = (pres[ZZ][ZZ]-(pres[XX][XX]+pres[YY][YY])*0.5)*box[ZZ][ZZ];
675 add_ebin(md->ebin,md->isurft,1,&tmp,bSum);
676 if (md->epc == epcPARRINELLORAHMAN || md->epc == epcMTTK)
677 {
678 tmp6[0] = state->boxv[XX][XX];
679 tmp6[1] = state->boxv[YY][YY];
680 tmp6[2] = state->boxv[ZZ][ZZ];
681 tmp6[3] = state->boxv[YY][XX];
682 tmp6[4] = state->boxv[ZZ][XX];
683 tmp6[5] = state->boxv[ZZ][YY];
684 add_ebin(md->ebin,md->ipc,md->bTricl ? 6 : 3,tmp6,bSum);
685 }
686 add_ebin(md->ebin,md->imu,3,mu_tot,bSum);
687 if (ekind && ekind->cosacc.cos_accel != 0)
688 {
689 vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
690 dens = (tmass*AMU)/(vol*NANO*NANO*NANO);
691 add_ebin(md->ebin,md->ivcos,1,&(ekind->cosacc.vcos),bSum);
692 /* 1/viscosity, unit 1/(kg m^-1 s^-1) */
693 tmp = 1/(ekind->cosacc.cos_accel/(ekind->cosacc.vcos*PICO)
694 *vol*sqr(box[ZZ][ZZ]*NANO/(2*M_PI)));
695 add_ebin(md->ebin,md->ivisc,1,&tmp,bSum);
696 }
697 if (md->nE > 1)
698 {
699 n=0;
700 for(i=0; (i<md->nEg); i++)
701 {
702 for(j=i; (j<md->nEg); j++)
703 {
704 gid=GID(i,j,md->nEg);
705 for(k=kk=0; (k<egNR); k++)
706 {
707 if (md->bEInd[k])
708 {
709 eee[kk++] = enerd->grpp.ener[k][gid];
710 }
711 }
712 add_ebin(md->ebin,md->igrp[n],md->nEc,eee,bSum);
713 n++;
714 }
715 }
716 }
717
718 if (ekind)
719 {
720 for(i=0; (i<md->nTC); i++)
721 {
722 md->tmp_r[i] = ekind->tcstat[i].T;
723 }
724 add_ebin(md->ebin,md->itemp,md->nTC,md->tmp_r,bSum);
725
726 bNoseHoover = (getenv("GMX_NOSEHOOVER_CHAINS") != NULL); /* whether to print Nose-Hoover chains */
727
728 if (md->etc == etcNOSEHOOVER)
729 {
730 if (bNoseHoover)
731 {
732 if (md->bNHC_trotter)
733 {
734 for(i=0; (i<md->nTC); i++)
735 {
736 for (j=0;j<md->nNHC;j++)
737 {
738 k = i*md->nNHC+j;
739 md->tmp_r[2*k] = state->nosehoover_xi[k];
740 md->tmp_r[2*k+1] = state->nosehoover_vxi[k];
741 }
742 }
743 add_ebin(md->ebin,md->itc,md->mde_n,md->tmp_r,bSum);
744
745 if (md->bMTTK) {
746 for(i=0; (i<md->nTCP); i++)
747 {
748 for (j=0;j<md->nNHC;j++)
749 {
750 k = i*md->nNHC+j;
751 md->tmp_r[2*k] = state->nhpres_xi[k];
752 md->tmp_r[2*k+1] = state->nhpres_vxi[k];
753 }
754 }
755 add_ebin(md->ebin,md->itcb,md->mdeb_n,md->tmp_r,bSum);
756 }
757
758 }
759 else
760 {
761 for(i=0; (i<md->nTC); i++)
762 {
763 md->tmp_r[2*i] = state->nosehoover_xi[i];
764 md->tmp_r[2*i+1] = state->nosehoover_vxi[i];
765 }
766 add_ebin(md->ebin,md->itc,md->mde_n,md->tmp_r,bSum);
767 }
768 }
769 }
770 else if (md->etc == etcBERENDSEN || md->etc == etcYES || md->etc == etcVRESCALE)
771 {
772 for(i=0; (i<md->nTC); i++)
773 {
774 md->tmp_r[i] = ekind->tcstat[i].lambda;
775 }
776 add_ebin(md->ebin,md->itc,md->nTC,md->tmp_r,bSum);
777 }
778 }
779
780 if (ekind && md->nU > 1)
781 {
782 for(i=0; (i<md->nU); i++)
783 {
784 copy_rvec(ekind->grpstat[i].u,md->tmp_v[i]);
785 }
786 add_ebin(md->ebin,md->iu,3*md->nU,md->tmp_v[0],bSum);
787 }
788
789 ebin_increase_count(md->ebin,bSum);
790
791 if (fp_dhdl)
792 {
793 fprintf(fp_dhdl,"%.4f %g",
794 time,
795 enerd->term[F_DVDL]+enerd->term[F_DKDL]+enerd->term[F_DHDL_CON]);
796 for(i=1; i<enerd->n_lambda; i++)
797 {
798 fprintf(fp_dhdl," %g",
799 enerd->enerpart_lambda[i]-enerd->enerpart_lambda[0]);
800 }
801 fprintf(fp_dhdl,"\n");
802 }
803 }
804
805 void upd_mdebin_step(t_mdebin *md)
806 {
807 ebin_increase_count(md->ebin,FALSE);
808 }
809
810 static void npr(FILE *log,int n,char c)
811 {
812 for(; (n>0); n--) fprintf(log,"%c",c);
813 }
814
815 static void pprint(FILE *log,const char *s,t_mdebin *md)
816 {
817 char CHAR='#';
818 int slen;
819 char buf1[22],buf2[22];
820
821 slen = strlen(s);
822 fprintf(log,"\t<====== ");
823 npr(log,slen,CHAR);
824 fprintf(log," ==>\n");
825 fprintf(log,"\t<==== %s ====>\n",s);
826 fprintf(log,"\t<== ");
827 npr(log,slen,CHAR);
828 fprintf(log," ======>\n\n");
829
830 fprintf(log,"\tStatistics over %s steps using %s frames\n",
831 gmx_step_str(md->ebin->nsteps_sim,buf1),
832 gmx_step_str(md->ebin->nsum_sim,buf2));
833 fprintf(log,"\n");
834 }
835
836 void print_ebin_header(FILE *log,gmx_large_int_t steps,double time,real lamb)
837 {
838 char buf[22];
839
840 fprintf(log," %12s %12s %12s\n"
841 " %12s %12.5f %12.5f\n\n",
842 "Step","Time","Lambda",gmx_step_str(steps,buf),time,lamb);
843 }
844
845 void print_ebin(ener_file_t fp_ene,bool bEne,bool bDR,bool bOR,
846 FILE *log,
847 gmx_large_int_t step,double time,
848 int mode,bool bCompact,
849 t_mdebin *md,t_fcdata *fcd,
850 gmx_groups_t *groups,t_grpopts *opts)
851 {
852 /*static char **grpnms=NULL;*/
853 char buf[246];
854 int i,j,n,ni,nj,ndr,nor;
855 int nr[enxNR];
856 real *block[enxNR];
857 t_enxframe fr;
858
859 switch (mode)
860 {
861 case eprNORMAL:
862 fr.t = time;
863 fr.step = step;
864 fr.nsteps = md->ebin->nsteps;
865 fr.nsum = md->ebin->nsum;
866 fr.nre = (bEne) ? md->ebin->nener : 0;
867 fr.ener = md->ebin->e;
868 fr.ndisre = bDR ? fcd->disres.npair : 0;
869 fr.disre_rm3tav = fcd->disres.rm3tav;
870 fr.disre_rt = fcd->disres.rt;
871 /* Optional additional blocks */
872 for(i=0; i<enxNR; i++)
873 {
874 nr[i] = 0;
875 }
876 if (fcd->orires.nr > 0 && bOR)
877 {
878 diagonalize_orires_tensors(&(fcd->orires));
879 nr[enxOR] = fcd->orires.nr;
880 block[enxOR] = fcd->orires.otav;
881 nr[enxORI] = (fcd->orires.oinsl != fcd->orires.otav) ?
882 fcd->orires.nr : 0;
883 block[enxORI] = fcd->orires.oinsl;
884 nr[enxORT] = fcd->orires.nex*12;
885 block[enxORT] = fcd->orires.eig;
886 }
887 fr.nblock = 0;
888 for(i=0; i<enxNR; i++)
889 {
890 if (nr[i] > 0)
891 {
892 fr.nblock = i + 1;
893 }
894 }
895 fr.nr = nr;
896 fr.block = block;
897 if (fr.nre || fr.ndisre || fr.nr[enxOR] || fr.nr[enxORI])
898 {
899 do_enx(fp_ene,&fr);
900 gmx_fio_check_file_position(enx_file_pointer(fp_ene));
901 if (fr.nre)
902 {
903 /* We have stored the sums, so reset the sum history */
904 reset_ebin_sums(md->ebin);
905 }
906 }
907 break;
908 case eprAVER:
909 if (log)
910 {
911 pprint(log,"A V E R A G E S",md);
912 }
913 break;
914 case eprRMS:
915 if (log)
916 {
917 pprint(log,"R M S - F L U C T U A T I O N S",md);
918 }
919 break;
920 default:
921 gmx_fatal(FARGS,"Invalid print mode (%d)",mode);
922 }
923
924 if (log)
925 {
926 for(i=0;i<opts->ngtc;i++)
927 {
928 if(opts->annealing[i]!=eannNO)
929 {
930 fprintf(log,"Current ref_t for group %s: %8.1f\n",
931 *(groups->grpname[groups->grps[egcTC].nm_ind[i]]),
932 opts->ref_t[i]);
933 }
934 }
935 if (mode==eprNORMAL && fcd->orires.nr>0)
936 {
937 print_orires_log(log,&(fcd->orires));
938 }
939 fprintf(log," Energies (%s)\n",unit_energy);
940 pr_ebin(log,md->ebin,md->ie,md->f_nre+md->nCrmsd,5,mode,TRUE);
941 fprintf(log,"\n");
942
943 if (!bCompact)
944 {
945 if (md->bDynBox)
946 {
947 pr_ebin(log,md->ebin,md->ib, md->bTricl ? NTRICLBOXS : NBOXS,5,
948 mode,TRUE);
949 fprintf(log,"\n");
950 }
951 if (md->bConstrVir)
952 {
953 fprintf(log," Constraint Virial (%s)\n",unit_energy);
954 pr_ebin(log,md->ebin,md->isvir,9,3,mode,FALSE);
955 fprintf(log,"\n");
956 fprintf(log," Force Virial (%s)\n",unit_energy);
957 pr_ebin(log,md->ebin,md->ifvir,9,3,mode,FALSE);
958 fprintf(log,"\n");
959 }
960 fprintf(log," Total Virial (%s)\n",unit_energy);
961 pr_ebin(log,md->ebin,md->ivir,9,3,mode,FALSE);
962 fprintf(log,"\n");
963 fprintf(log," Pressure (%s)\n",unit_pres_bar);
964 pr_ebin(log,md->ebin,md->ipres,9,3,mode,FALSE);
965 fprintf(log,"\n");
966 fprintf(log," Total Dipole (%s)\n",unit_dipole_D);
967 pr_ebin(log,md->ebin,md->imu,3,3,mode,FALSE);
968 fprintf(log,"\n");
969
970 if (md->nE > 1)
971 {
972 if (md->print_grpnms==NULL)
973 {
974 snew(md->print_grpnms,md->nE);
975 n=0;
976 for(i=0; (i<md->nEg); i++)
977 {
978 ni=groups->grps[egcENER].nm_ind[i];
979 for(j=i; (j<md->nEg); j++)
980 {
981 nj=groups->grps[egcENER].nm_ind[j];
982 sprintf(buf,"%s-%s",*(groups->grpname[ni]),
983 *(groups->grpname[nj]));
984 md->print_grpnms[n++]=strdup(buf);
985 }
986 }
987 }
988 sprintf(buf,"Epot (%s)",unit_energy);
989 fprintf(log,"%15s ",buf);
990 for(i=0; (i<egNR); i++)
991 {
992 if (md->bEInd[i])
993 {
994 fprintf(log,"%12s ",egrp_nm[i]);
995 }
996 }
997 fprintf(log,"\n");
998 for(i=0; (i<md->nE); i++)
999 {
1000 fprintf(log,"%15s",md->print_grpnms[i]);
1001 pr_ebin(log,md->ebin,md->igrp[i],md->nEc,md->nEc,mode,
1002 FALSE);
1003 }
1004 fprintf(log,"\n");
1005 }
1006 if (md->nTC > 1)
1007 {
1008 pr_ebin(log,md->ebin,md->itemp,md->nTC,4,mode,TRUE);
1009 fprintf(log,"\n");
1010 }
1011 if (md->nU > 1)
1012 {
1013 fprintf(log,"%15s %12s %12s %12s\n",
1014 "Group","Ux","Uy","Uz");
1015 for(i=0; (i<md->nU); i++)
1016 {
1017 ni=groups->grps[egcACC].nm_ind[i];
1018 fprintf(log,"%15s",*groups->grpname[ni]);
1019 pr_ebin(log,md->ebin,md->iu+3*i,3,3,mode,FALSE);
1020 }
1021 fprintf(log,"\n");
1022 }
1023 }
1024
1025 }
1026 }
1027
1028 void
1029 init_energyhistory(energyhistory_t * enerhist)
1030 {
1031 enerhist->nener = 0;
1032
1033 enerhist->ener_ave = NULL;
1034 enerhist->ener_sum = NULL;
1035 enerhist->ener_sum_sim = NULL;
1036
1037 enerhist->nsteps = 0;
1038 enerhist->nsum = 0;
1039 enerhist->nsteps_sim = 0;
1040 enerhist->nsum_sim = 0;
1041 }
1042
1043 void
1044 update_energyhistory(energyhistory_t * enerhist,t_mdebin * mdebin)
1045 {
1046 int i;
1047
1048 enerhist->nsteps = mdebin->ebin->nsteps;
1049 enerhist->nsum = mdebin->ebin->nsum;
1050 enerhist->nsteps_sim = mdebin->ebin->nsteps_sim;
1051 enerhist->nsum_sim = mdebin->ebin->nsum_sim;
1052 enerhist->nener = mdebin->ebin->nener;
1053
1054 if (mdebin->ebin->nsum > 0)
1055 {
1056 /* Check if we need to allocate first */
1057 if(enerhist->ener_ave == NULL)
1058 {
1059 snew(enerhist->ener_ave,enerhist->nener);
1060 snew(enerhist->ener_sum,enerhist->nener);
1061 }
1062
1063 for(i=0;i<enerhist->nener;i++)
1064 {
1065 enerhist->ener_ave[i] = mdebin->ebin->e[i].eav;
1066 enerhist->ener_sum[i] = mdebin->ebin->e[i].esum;
1067 }
1068 }
1069
1070 if (mdebin->ebin->nsum_sim > 0)
1071 {
1072 /* Check if we need to allocate first */
1073 if(enerhist->ener_sum_sim == NULL)
1074 {
1075 snew(enerhist->ener_sum_sim,enerhist->nener);
1076 }
1077
1078 for(i=0;i<enerhist->nener;i++)
1079 {
1080 enerhist->ener_sum_sim[i] = mdebin->ebin->e_sim[i].esum;
1081 }
1082 }
1083 }
1084
1085 void
1086 restore_energyhistory_from_state(t_mdebin * mdebin,energyhistory_t * enerhist)
1087 {
1088 int i;
1089
1090 if ((enerhist->nsum > 0 || enerhist->nsum_sim > 0) &&
1091 mdebin->ebin->nener != enerhist->nener)
1092 {
1093 gmx_fatal(FARGS,"Mismatch between number of energies in run input (%d) and checkpoint file (%d).",
1094 mdebin->ebin->nener,enerhist->nener);
1095 }
1096
1097 mdebin->ebin->nsteps = enerhist->nsteps;
1098 mdebin->ebin->nsum = enerhist->nsum;
1099 mdebin->ebin->nsteps_sim = enerhist->nsteps_sim;
1100 mdebin->ebin->nsum_sim = enerhist->nsum_sim;
1101
1102 for(i=0; i<mdebin->ebin->nener; i++)
1103 {
1104 mdebin->ebin->e[i].eav =
1105 (enerhist->nsum > 0 ? enerhist->ener_ave[i] : 0);
1106 mdebin->ebin->e[i].esum =
1107 (enerhist->nsum > 0 ? enerhist->ener_sum[i] : 0);
1108 mdebin->ebin->e_sim[i].esum =
1109 (enerhist->nsum_sim > 0 ? enerhist->ener_sum_sim[i] : 0);
1110 }
1111 }
Atempt to write QMMM interface for Gamess-US