?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
added Verlet scheme and NxN non-bonded functionality
[gromacs.git] / src / gmxlib / typedefs.c
blob7bec9cd9994809c6b5843a666827d417094330bc
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 * GROningen Mixture of Alchemy and Childrens' Stories
35 */
36 /* This file is completely threadsafe - keep it that way! */
37 #ifdef HAVE_CONFIG_H
38 #include <config.h>
39 #endif
40
41 #include "smalloc.h"
42 #include "symtab.h"
43 #include "vec.h"
44 #include "pbc.h"
45 #include <string.h>
46
47 #ifdef GMX_THREAD_MPI
48 #include "thread_mpi.h"
49 #endif
50
51 /* The source code in this file should be thread-safe.
52 Please keep it that way. */
53
54
55
56 static gmx_bool bOverAllocDD=FALSE;
57 #ifdef GMX_THREAD_MPI
58 static tMPI_Thread_mutex_t over_alloc_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
59 #endif
60
61
62 void set_over_alloc_dd(gmx_bool set)
63 {
64 #ifdef GMX_THREAD_MPI
65 tMPI_Thread_mutex_lock(&over_alloc_mutex);
66 /* we just make sure that we don't set this at the same time.
67 We don't worry too much about reading this rarely-set variable */
68 #endif
69 bOverAllocDD = set;
70 #ifdef GMX_THREAD_MPI
71 tMPI_Thread_mutex_unlock(&over_alloc_mutex);
72 #endif
73 }
74
75 int over_alloc_dd(int n)
76 {
77 if (bOverAllocDD)
78 return OVER_ALLOC_FAC*n + 100;
79 else
80 return n;
81 }
82
83 int gmx_large_int_to_int(gmx_large_int_t step,const char *warn)
84 {
85 int i;
86
87 i = (int)step;
88
89 if (warn != NULL && (step < INT_MIN || step > INT_MAX)) {
90 fprintf(stderr,"\nWARNING during %s:\n",warn);
91 fprintf(stderr,"step value ");
92 fprintf(stderr,gmx_large_int_pfmt,step);
93 fprintf(stderr," does not fit in int, converted to %d\n\n",i);
94 }
95
96 return i;
97 }
98
99 char *gmx_step_str(gmx_large_int_t i,char *buf)
100 {
101 sprintf(buf,gmx_large_int_pfmt,i);
102
103 return buf;
104 }
105
106 void init_block(t_block *block)
107 {
108 int i;
109
110 block->nr = 0;
111 block->nalloc_index = 1;
112 snew(block->index,block->nalloc_index);
113 block->index[0] = 0;
114 }
115
116 void init_blocka(t_blocka *block)
117 {
118 int i;
119
120 block->nr = 0;
121 block->nra = 0;
122 block->nalloc_index = 1;
123 snew(block->index,block->nalloc_index);
124 block->index[0] = 0;
125 block->nalloc_a = 0;
126 block->a = NULL;
127 }
128
129 void init_atom(t_atoms *at)
130 {
131 int i;
132
133 at->nr = 0;
134 at->nres = 0;
135 at->atom = NULL;
136 at->resinfo = NULL;
137 at->atomname = NULL;
138 at->atomtype = NULL;
139 at->atomtypeB= NULL;
140 at->pdbinfo = NULL;
141 }
142
143 void init_atomtypes(t_atomtypes *at)
144 {
145 at->nr = 0;
146 at->radius = NULL;
147 at->vol = NULL;
148 at->atomnumber = NULL;
149 at->gb_radius = NULL;
150 at->S_hct = NULL;
151 }
152
153 void init_groups(gmx_groups_t *groups)
154 {
155 int g;
156
157 groups->ngrpname = 0;
158 groups->grpname = NULL;
159 for(g=0; (g<egcNR); g++) {
160 groups->grps[g].nm_ind = NULL;
161 groups->ngrpnr[g] = 0;
162 groups->grpnr[g] = NULL;
163 }
164
165 }
166
167 void init_mtop(gmx_mtop_t *mtop)
168 {
169 mtop->name = NULL;
170 mtop->nmoltype = 0;
171 mtop->moltype = NULL;
172 mtop->nmolblock = 0;
173 mtop->molblock = NULL;
174 mtop->maxres_renum = 0;
175 mtop->maxresnr = -1;
176 init_groups(&mtop->groups);
177 init_block(&mtop->mols);
178 open_symtab(&mtop->symtab);
179 }
180
181 void init_top (t_topology *top)
182 {
183 int i;
184
185 top->name = NULL;
186 init_atom (&(top->atoms));
187 init_atomtypes(&(top->atomtypes));
188 init_block(&top->cgs);
189 init_block(&top->mols);
190 init_blocka(&top->excls);
191 open_symtab(&top->symtab);
192 }
193
194 void init_inputrec(t_inputrec *ir)
195 {
196 memset(ir,0,(size_t)sizeof(*ir));
197 snew(ir->fepvals,1);
198 snew(ir->expandedvals,1);
199 snew(ir->simtempvals,1);
200 }
201
202 void stupid_fill_block(t_block *grp,int natom,gmx_bool bOneIndexGroup)
203 {
204 int i;
205
206 if (bOneIndexGroup) {
207 grp->nalloc_index = 2;
208 snew(grp->index,grp->nalloc_index);
209 grp->index[0]=0;
210 grp->index[1]=natom;
211 grp->nr=1;
212 }
213 else {
214 grp->nalloc_index = natom+1;
215 snew(grp->index,grp->nalloc_index);
216 snew(grp->index,natom+1);
217 for(i=0; (i<=natom); i++)
218 grp->index[i]=i;
219 grp->nr=natom;
220 }
221 }
222
223 void stupid_fill_blocka(t_blocka *grp,int natom)
224 {
225 int i;
226
227 grp->nalloc_a = natom;
228 snew(grp->a,grp->nalloc_a);
229 for(i=0; (i<natom); i++)
230 grp->a[i]=i;
231 grp->nra=natom;
232
233 grp->nalloc_index = natom + 1;
234 snew(grp->index,grp->nalloc_index);
235 for(i=0; (i<=natom); i++)
236 grp->index[i]=i;
237 grp->nr=natom;
238 }
239
240 void copy_blocka(const t_blocka *src,t_blocka *dest)
241 {
242 int i;
243
244 dest->nr = src->nr;
245 dest->nalloc_index = dest->nr + 1;
246 snew(dest->index,dest->nalloc_index);
247 for(i=0; i<dest->nr+1; i++) {
248 dest->index[i] = src->index[i];
249 }
250 dest->nra = src->nra;
251 dest->nalloc_a = dest->nra + 1;
252 snew(dest->a,dest->nalloc_a);
253 for(i=0; i<dest->nra+1; i++) {
254 dest->a[i] = src->a[i];
255 }
256 }
257
258 void done_block(t_block *block)
259 {
260 block->nr = 0;
261 sfree(block->index);
262 block->nalloc_index = 0;
263 }
264
265 void done_blocka(t_blocka *block)
266 {
267 block->nr = 0;
268 block->nra = 0;
269 sfree(block->index);
270 if (block->a)
271 sfree(block->a);
272 block->nalloc_index = 0;
273 block->nalloc_a = 0;
274 }
275
276 void done_atom (t_atoms *at)
277 {
278 at->nr = 0;
279 at->nres = 0;
280 sfree(at->atom);
281 sfree(at->resinfo);
282 sfree(at->atomname);
283 sfree(at->atomtype);
284 sfree(at->atomtypeB);
285 if (at->pdbinfo)
286 sfree(at->pdbinfo);
287 }
288
289 void done_atomtypes(t_atomtypes *atype)
290 {
291 atype->nr = 0;
292 sfree(atype->radius);
293 sfree(atype->vol);
294 sfree(atype->surftens);
295 sfree(atype->atomnumber);
296 sfree(atype->gb_radius);
297 sfree(atype->S_hct);
298 }
299
300 void done_moltype(gmx_moltype_t *molt)
301 {
302 int f;
303
304 done_atom(&molt->atoms);
305 done_block(&molt->cgs);
306 done_blocka(&molt->excls);
307
308 for(f=0; f<F_NRE; f++) {
309 sfree(molt->ilist[f].iatoms);
310 molt->ilist[f].nalloc = 0;
311 }
312 }
313
314 void done_molblock(gmx_molblock_t *molb)
315 {
316 if (molb->nposres_xA > 0) {
317 molb->nposres_xA = 0;
318 free(molb->posres_xA);
319 }
320 if (molb->nposres_xB > 0) {
321 molb->nposres_xB = 0;
322 free(molb->posres_xB);
323 }
324 }
325
326 void done_mtop(gmx_mtop_t *mtop,gmx_bool bDoneSymtab)
327 {
328 int i;
329
330 if (bDoneSymtab) {
331 done_symtab(&mtop->symtab);
332 }
333
334 sfree(mtop->ffparams.functype);
335 sfree(mtop->ffparams.iparams);
336
337 for(i=0; i<mtop->nmoltype; i++) {
338 done_moltype(&mtop->moltype[i]);
339 }
340 sfree(mtop->moltype);
341 for(i=0; i<mtop->nmolblock; i++) {
342 done_molblock(&mtop->molblock[i]);
343 }
344 sfree(mtop->molblock);
345 done_block(&mtop->mols);
346 }
347
348 void done_top(t_topology *top)
349 {
350 int f;
351
352 sfree(top->idef.functype);
353 sfree(top->idef.iparams);
354 for (f = 0; f < F_NRE; ++f)
355 {
356 sfree(top->idef.il[f].iatoms);
357 top->idef.il[f].iatoms = NULL;
358 top->idef.il[f].nalloc = 0;
359 }
360
361 done_atom (&(top->atoms));
362
363 /* For GB */
364 done_atomtypes(&(top->atomtypes));
365
366 done_symtab(&(top->symtab));
367 done_block(&(top->cgs));
368 done_block(&(top->mols));
369 done_blocka(&(top->excls));
370 }
371
372 static void done_pullgrp(t_pullgrp *pgrp)
373 {
374 sfree(pgrp->ind);
375 sfree(pgrp->ind_loc);
376 sfree(pgrp->weight);
377 sfree(pgrp->weight_loc);
378 }
379
380 static void done_pull(t_pull *pull)
381 {
382 int i;
383
384 for(i=0; i<pull->ngrp+1; i++) {
385 done_pullgrp(pull->grp);
386 done_pullgrp(pull->dyna);
387 }
388 }
389
390 void done_inputrec(t_inputrec *ir)
391 {
392 int m;
393
394 for(m=0; (m<DIM); m++) {
395 if (ir->ex[m].a) sfree(ir->ex[m].a);
396 if (ir->ex[m].phi) sfree(ir->ex[m].phi);
397 if (ir->et[m].a) sfree(ir->et[m].a);
398 if (ir->et[m].phi) sfree(ir->et[m].phi);
399 }
400
401 sfree(ir->opts.nrdf);
402 sfree(ir->opts.ref_t);
403 sfree(ir->opts.annealing);
404 sfree(ir->opts.anneal_npoints);
405 sfree(ir->opts.anneal_time);
406 sfree(ir->opts.anneal_temp);
407 sfree(ir->opts.tau_t);
408 sfree(ir->opts.acc);
409 sfree(ir->opts.nFreeze);
410 sfree(ir->opts.QMmethod);
411 sfree(ir->opts.QMbasis);
412 sfree(ir->opts.QMcharge);
413 sfree(ir->opts.QMmult);
414 sfree(ir->opts.bSH);
415 sfree(ir->opts.CASorbitals);
416 sfree(ir->opts.CASelectrons);
417 sfree(ir->opts.SAon);
418 sfree(ir->opts.SAoff);
419 sfree(ir->opts.SAsteps);
420 sfree(ir->opts.bOPT);
421 sfree(ir->opts.bTS);
422
423 if (ir->pull) {
424 done_pull(ir->pull);
425 sfree(ir->pull);
426 }
427 }
428
429 static void zero_ekinstate(ekinstate_t *eks)
430 {
431 eks->ekin_n = 0;
432 eks->ekinh = NULL;
433 eks->ekinf = NULL;
434 eks->ekinh_old = NULL;
435 eks->ekinscalef_nhc = NULL;
436 eks->ekinscaleh_nhc = NULL;
437 eks->vscale_nhc = NULL;
438 eks->dekindl = 0;
439 eks->mvcos = 0;
440 }
441
442 void init_energyhistory(energyhistory_t * enerhist)
443 {
444 enerhist->nener = 0;
445
446 enerhist->ener_ave = NULL;
447 enerhist->ener_sum = NULL;
448 enerhist->ener_sum_sim = NULL;
449 enerhist->dht = NULL;
450
451 enerhist->nsteps = 0;
452 enerhist->nsum = 0;
453 enerhist->nsteps_sim = 0;
454 enerhist->nsum_sim = 0;
455
456 enerhist->dht = NULL;
457 }
458
459 static void done_delta_h_history(delta_h_history_t *dht)
460 {
461 int i;
462
463 for(i=0; i<dht->nndh; i++)
464 {
465 sfree(dht->dh[i]);
466 }
467 sfree(dht->dh);
468 sfree(dht->ndh);
469 }
470
471 void done_energyhistory(energyhistory_t * enerhist)
472 {
473 sfree(enerhist->ener_ave);
474 sfree(enerhist->ener_sum);
475 sfree(enerhist->ener_sum_sim);
476
477 if (enerhist->dht != NULL)
478 {
479 done_delta_h_history(enerhist->dht);
480 sfree(enerhist->dht);
481 }
482 }
483
484 void init_gtc_state(t_state *state, int ngtc, int nnhpres, int nhchainlength)
485 {
486 int i,j;
487
488 state->ngtc = ngtc;
489 state->nnhpres = nnhpres;
490 state->nhchainlength = nhchainlength;
491 if (state->ngtc > 0)
492 {
493 snew(state->nosehoover_xi,state->nhchainlength*state->ngtc);
494 snew(state->nosehoover_vxi,state->nhchainlength*state->ngtc);
495 snew(state->therm_integral,state->ngtc);
496 for(i=0; i<state->ngtc; i++)
497 {
498 for (j=0;j<state->nhchainlength;j++)
499 {
500 state->nosehoover_xi[i*state->nhchainlength + j] = 0.0;
501 state->nosehoover_vxi[i*state->nhchainlength + j] = 0.0;
502 }
503 }
504 for(i=0; i<state->ngtc; i++) {
505 state->therm_integral[i] = 0.0;
506 }
507 }
508 else
509 {
510 state->nosehoover_xi = NULL;
511 state->nosehoover_vxi = NULL;
512 state->therm_integral = NULL;
513 }
514
515 if (state->nnhpres > 0)
516 {
517 snew(state->nhpres_xi,state->nhchainlength*nnhpres);
518 snew(state->nhpres_vxi,state->nhchainlength*nnhpres);
519 for(i=0; i<nnhpres; i++)
520 {
521 for (j=0;j<state->nhchainlength;j++)
522 {
523 state->nhpres_xi[i*nhchainlength + j] = 0.0;
524 state->nhpres_vxi[i*nhchainlength + j] = 0.0;
525 }
526 }
527 }
528 else
529 {
530 state->nhpres_xi = NULL;
531 state->nhpres_vxi = NULL;
532 }
533 }
534
535
536 void init_state(t_state *state, int natoms, int ngtc, int nnhpres, int nhchainlength, int nlambda)
537 {
538 int i;
539
540 state->natoms = natoms;
541 state->nrng = 0;
542 state->flags = 0;
543 state->lambda = 0;
544 snew(state->lambda,efptNR);
545 for (i=0;i<efptNR;i++)
546 {
547 state->lambda[i] = 0;
548 }
549 state->veta = 0;
550 clear_mat(state->box);
551 clear_mat(state->box_rel);
552 clear_mat(state->boxv);
553 clear_mat(state->pres_prev);
554 clear_mat(state->svir_prev);
555 clear_mat(state->fvir_prev);
556 init_gtc_state(state,ngtc,nnhpres,nhchainlength);
557 state->nalloc = state->natoms;
558 if (state->nalloc > 0) {
559 snew(state->x,state->nalloc);
560 snew(state->v,state->nalloc);
561 } else {
562 state->x = NULL;
563 state->v = NULL;
564 }
565 state->sd_X = NULL;
566 state->cg_p = NULL;
567
568 zero_ekinstate(&state->ekinstate);
569
570 init_energyhistory(&state->enerhist);
571
572 init_df_history(&state->dfhist,nlambda,0);
573
574 state->ddp_count = 0;
575 state->ddp_count_cg_gl = 0;
576 state->cg_gl = NULL;
577 state->cg_gl_nalloc = 0;
578 }
579
580 void done_state(t_state *state)
581 {
582 if (state->nosehoover_xi) sfree(state->nosehoover_xi);
583 if (state->x) sfree(state->x);
584 if (state->v) sfree(state->v);
585 if (state->sd_X) sfree(state->sd_X);
586 if (state->cg_p) sfree(state->cg_p);
587 state->nalloc = 0;
588 if (state->cg_gl) sfree(state->cg_gl);
589 state->cg_gl_nalloc = 0;
590 }
591
592 static void do_box_rel(t_inputrec *ir,matrix box_rel,matrix b,gmx_bool bInit)
593 {
594 int d,d2;
595
596 for(d=YY; d<=ZZ; d++) {
597 for(d2=XX; d2<=(ir->epct==epctSEMIISOTROPIC ? YY : ZZ); d2++) {
598 /* We need to check if this box component is deformed
599 * or if deformation of another component might cause
600 * changes in this component due to box corrections.
601 */
602 if (ir->deform[d][d2] == 0 &&
603 !(d == ZZ && d2 == XX && ir->deform[d][YY] != 0 &&
604 (b[YY][d2] != 0 || ir->deform[YY][d2] != 0))) {
605 if (bInit) {
606 box_rel[d][d2] = b[d][d2]/b[XX][XX];
607 } else {
608 b[d][d2] = b[XX][XX]*box_rel[d][d2];
609 }
610 }
611 }
612 }
613 }
614
615 void set_box_rel(t_inputrec *ir,t_state *state)
616 {
617 /* Make sure the box obeys the restrictions before we fix the ratios */
618 correct_box(NULL,0,state->box,NULL);
619
620 clear_mat(state->box_rel);
621
622 if (PRESERVE_SHAPE(*ir))
623 do_box_rel(ir,state->box_rel,state->box,TRUE);
624 }
625
626 void preserve_box_shape(t_inputrec *ir,matrix box_rel,matrix b)
627 {
628 if (PRESERVE_SHAPE(*ir))
629 do_box_rel(ir,box_rel,b,FALSE);
630 }
631
632 void add_t_atoms(t_atoms *atoms,int natom_extra,int nres_extra)
633 {
634 int i;
635
636 if (natom_extra > 0)
637 {
638 srenew(atoms->atomname,atoms->nr+natom_extra);
639 srenew(atoms->atom,atoms->nr+natom_extra);
640 if (NULL != atoms->pdbinfo)
641 srenew(atoms->pdbinfo,atoms->nr+natom_extra);
642 if (NULL != atoms->atomtype)
643 srenew(atoms->atomtype,atoms->nr+natom_extra);
644 if (NULL != atoms->atomtypeB)
645 srenew(atoms->atomtypeB,atoms->nr+natom_extra);
646 for(i=atoms->nr; (i<atoms->nr+natom_extra); i++) {
647 atoms->atomname[i] = NULL;
648 memset(&atoms->atom[i],0,sizeof(atoms->atom[i]));
649 if (NULL != atoms->pdbinfo)
650 memset(&atoms->pdbinfo[i],0,sizeof(atoms->pdbinfo[i]));
651 if (NULL != atoms->atomtype)
652 atoms->atomtype[i] = NULL;
653 if (NULL != atoms->atomtypeB)
654 atoms->atomtypeB[i] = NULL;
655 }
656 atoms->nr += natom_extra;
657 }
658 if (nres_extra > 0)
659 {
660 srenew(atoms->resinfo,atoms->nres+nres_extra);
661 for(i=atoms->nres; (i<atoms->nres+nres_extra); i++) {
662 memset(&atoms->resinfo[i],0,sizeof(atoms->resinfo[i]));
663 }
664 atoms->nres += nres_extra;
665 }
666 }
667
668 void init_t_atoms(t_atoms *atoms, int natoms, gmx_bool bPdbinfo)
669 {
670 atoms->nr=natoms;
671 atoms->nres=0;
672 snew(atoms->atomname,natoms);
673 atoms->atomtype=NULL;
674 atoms->atomtypeB=NULL;
675 snew(atoms->resinfo,natoms);
676 snew(atoms->atom,natoms);
677 if (bPdbinfo)
678 snew(atoms->pdbinfo,natoms);
679 else
680 atoms->pdbinfo=NULL;
681 }
682
683 t_atoms *copy_t_atoms(t_atoms *src)
684 {
685 t_atoms *dst;
686 int i;
687
688 snew(dst,1);
689 init_t_atoms(dst,src->nr,(NULL != src->pdbinfo));
690 dst->nr = src->nr;
691 if (NULL != src->atomname)
692 snew(dst->atomname,src->nr);
693 if (NULL != src->atomtype)
694 snew(dst->atomtype,src->nr);
695 if (NULL != src->atomtypeB)
696 snew(dst->atomtypeB,src->nr);
697 for(i=0; (i<src->nr); i++) {
698 dst->atom[i] = src->atom[i];
699 if (NULL != src->pdbinfo)
700 dst->pdbinfo[i] = src->pdbinfo[i];
701 if (NULL != src->atomname)
702 dst->atomname[i] = src->atomname[i];
703 if (NULL != src->atomtype)
704 dst->atomtype[i] = src->atomtype[i];
705 if (NULL != src->atomtypeB)
706 dst->atomtypeB[i] = src->atomtypeB[i];
707 }
708 dst->nres = src->nres;
709 for(i=0; (i<src->nres); i++) {
710 dst->resinfo[i] = src->resinfo[i];
711 }
712 return dst;
713 }
714
715 void t_atoms_set_resinfo(t_atoms *atoms,int atom_ind,t_symtab *symtab,
716 const char *resname,int resnr,unsigned char ic,
717 int chainnum, char chainid)
718 {
719 t_resinfo *ri;
720
721 ri = &atoms->resinfo[atoms->atom[atom_ind].resind];
722 ri->name = put_symtab(symtab,resname);
723 ri->rtp = NULL;
724 ri->nr = resnr;
725 ri->ic = ic;
726 ri->chainnum = chainnum;
727 ri->chainid = chainid;
728 }
729
730 void free_t_atoms(t_atoms *atoms,gmx_bool bFreeNames)
731 {
732 int i;
733
734 if (bFreeNames) {
735 for(i=0; i<atoms->nr; i++) {
736 sfree(*atoms->atomname[i]);
737 *atoms->atomname[i]=NULL;
738 }
739 for(i=0; i<atoms->nres; i++) {
740 sfree(*atoms->resinfo[i].name);
741 *atoms->resinfo[i].name=NULL;
742 }
743 }
744 sfree(atoms->atomname);
745 /* Do we need to free atomtype and atomtypeB as well ? */
746 sfree(atoms->resinfo);
747 sfree(atoms->atom);
748 if (atoms->pdbinfo)
749 sfree(atoms->pdbinfo);
750 atoms->nr=0;
751 atoms->nres=0;
752 }
753
754 real max_cutoff(real cutoff1,real cutoff2)
755 {
756 if (cutoff1 == 0 || cutoff2 == 0)
757 {
758 return 0;
759 }
760 else
761 {
762 return max(cutoff1,cutoff2);
763 }
764 }
765
766 extern void init_df_history(df_history_t *dfhist, int nlambda, real wl_delta)
767 {
768 int i;
769
770 dfhist->bEquil = 0;
771 dfhist->nlambda = nlambda;
772 dfhist->wl_delta = wl_delta;
773 snew(dfhist->sum_weights,dfhist->nlambda);
774 snew(dfhist->sum_dg,dfhist->nlambda);
775 snew(dfhist->sum_minvar,dfhist->nlambda);
776 snew(dfhist->sum_variance,dfhist->nlambda);
777 snew(dfhist->n_at_lam,dfhist->nlambda);
778 snew(dfhist->wl_histo,dfhist->nlambda);
779
780 /* allocate transition matrices here */
781 snew(dfhist->Tij,dfhist->nlambda);
782 snew(dfhist->Tij_empirical,dfhist->nlambda);
783
784 for (i=0;i<dfhist->nlambda;i++) {
785 snew(dfhist->Tij[i],dfhist->nlambda);
786 snew(dfhist->Tij_empirical[i],dfhist->nlambda);
787 }
788
789 snew(dfhist->accum_p,dfhist->nlambda);
790 snew(dfhist->accum_m,dfhist->nlambda);
791 snew(dfhist->accum_p2,dfhist->nlambda);
792 snew(dfhist->accum_m2,dfhist->nlambda);
793
794 for (i=0;i<dfhist->nlambda;i++) {
795 snew((dfhist->accum_p)[i],dfhist->nlambda);
796 snew((dfhist->accum_m)[i],dfhist->nlambda);
797 snew((dfhist->accum_p2)[i],dfhist->nlambda);
798 snew((dfhist->accum_m2)[i],dfhist->nlambda);
799 }
800 }
801
802 extern void copy_df_history(df_history_t *df_dest, df_history_t *df_source)
803 {
804 int i,j;
805
806 init_df_history(df_dest,df_source->nlambda,df_source->wl_delta);
807 df_dest->nlambda = df_source->nlambda;
808 df_dest->bEquil = df_source->bEquil;
809 for (i=0;i<df_dest->nlambda;i++)
810 {
811 df_dest->sum_weights[i] = df_source->sum_weights[i];
812 df_dest->sum_dg[i] = df_source->sum_dg[i];
813 df_dest->sum_minvar[i] = df_source->sum_minvar[i];
814 df_dest->sum_variance[i] = df_source->sum_variance[i];
815 df_dest->n_at_lam[i] = df_source->n_at_lam[i];
816 df_dest->wl_histo[i] = df_source->wl_histo[i];
817 df_dest->accum_p[i] = df_source->accum_p[i];
818 df_dest->accum_m[i] = df_source->accum_m[i];
819 df_dest->accum_p2[i] = df_source->accum_p2[i];
820 df_dest->accum_m2[i] = df_source->accum_m2[i];
821 }
822
823 for (i=0;i<df_dest->nlambda;i++)
824 {
825 for (j=0;j<df_dest->nlambda;j++)
826 {
827 df_dest->Tij[i][j] = df_source->Tij[i][j];
828 df_dest->Tij_empirical[i][j] = df_source->Tij_empirical[i][j];
829 }
830 }
831 }
The ultimate molecular dynamics simulation package