Project

General

Profile

typedefs.c

edited in order to initialize arrays - Christopher Mirabzadeh, 11/02/2015 11:14 PM

 
1
/*
2
 * This file is part of the GROMACS molecular simulation package.
3
 *
4
 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5
 * Copyright (c) 2001-2004, The GROMACS development team.
6
 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7
 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8
 * and including many others, as listed in the AUTHORS file in the
9
 * top-level source directory and at http://www.gromacs.org.
10
 *
11
 * GROMACS is free software; you can redistribute it and/or
12
 * modify it under the terms of the GNU Lesser General Public License
13
 * as published by the Free Software Foundation; either version 2.1
14
 * of the License, or (at your option) any later version.
15
 *
16
 * GROMACS is distributed in the hope that it will be useful,
17
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
19
 * Lesser General Public License for more details.
20
 *
21
 * You should have received a copy of the GNU Lesser General Public
22
 * License along with GROMACS; if not, see
23
 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
25
 *
26
 * If you want to redistribute modifications to GROMACS, please
27
 * consider that scientific software is very special. Version
28
 * control is crucial - bugs must be traceable. We will be happy to
29
 * consider code for inclusion in the official distribution, but
30
 * derived work must not be called official GROMACS. Details are found
31
 * in the README & COPYING files - if they are missing, get the
32
 * official version at http://www.gromacs.org.
33
 *
34
 * To help us fund GROMACS development, we humbly ask that you cite
35
 * the research papers on the package. Check out http://www.gromacs.org.
36
 */
37
/* This file is completely threadsafe - keep it that way! */
38
#ifdef HAVE_CONFIG_H
39
#include <config.h>
40
#endif
41

    
42
#include "thread_mpi/threads.h"
43

    
44
#include "gromacs/utility/smalloc.h"
45
#include "symtab.h"
46
#include "vec.h"
47
#include "pbc.h"
48
#include "macros.h"
49
#include <string.h>
50
#include "gromacs/random/random.h"
51

    
52
/* The source code in this file should be thread-safe.
53
      Please keep it that way. */
54

    
55

    
56

    
57
static gmx_bool            bOverAllocDD     = FALSE;
58
static tMPI_Thread_mutex_t over_alloc_mutex = TMPI_THREAD_MUTEX_INITIALIZER;
59

    
60

    
61
void set_over_alloc_dd(gmx_bool set)
62
{
63
    tMPI_Thread_mutex_lock(&over_alloc_mutex);
64
    /* we just make sure that we don't set this at the same time.
65
       We don't worry too much about reading this rarely-set variable */
66
    bOverAllocDD = set;
67
    tMPI_Thread_mutex_unlock(&over_alloc_mutex);
68
}
69

    
70
int over_alloc_dd(int n)
71
{
72
    if (bOverAllocDD)
73
    {
74
        return OVER_ALLOC_FAC*n + 100;
75
    }
76
    else
77
    {
78
        return n;
79
    }
80
}
81

    
82
int gmx_int64_to_int(gmx_int64_t step, const char *warn)
83
{
84
    int i;
85

    
86
    i = (int)step;
87

    
88
    if (warn != NULL && (step < INT_MIN || step > INT_MAX))
89
    {
90
        fprintf(stderr, "\nWARNING during %s:\n", warn);
91
        fprintf(stderr, "step value ");
92
        fprintf(stderr, "%"GMX_PRId64, 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_int64_t i, char *buf)
100
{
101
    sprintf(buf, "%"GMX_PRId64, 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
    {
161
        groups->grps[g].nm_ind = NULL;
162
        groups->ngrpnr[g]      = 0;
163
        groups->grpnr[g]       = NULL;
164
    }
165

    
166
}
167

    
168
void init_mtop(gmx_mtop_t *mtop)
169
{
170
    mtop->name         = NULL;
171
    mtop->nmoltype     = 0;
172
    mtop->moltype      = NULL;
173
    mtop->nmolblock    = 0;
174
    mtop->molblock     = NULL;
175
    mtop->maxres_renum = 0;
176
    mtop->maxresnr     = -1;
177
    init_groups(&mtop->groups);
178
    init_block(&mtop->mols);
179
    open_symtab(&mtop->symtab);
180
}
181

    
182
void init_top(t_topology *top)
183
{
184
    int i;
185

    
186
    top->name = NULL;
187
    init_atom (&(top->atoms));
188
    init_atomtypes(&(top->atomtypes));
189
    init_block(&top->cgs);
190
    init_block(&top->mols);
191
    init_blocka(&top->excls);
192
    open_symtab(&top->symtab);
193
}
194

    
195
void init_inputrec(t_inputrec *ir)
196
{
197
    memset(ir, 0, (size_t)sizeof(*ir));
198
    snew(ir->fepvals, 1);
199
    snew(ir->expandedvals, 1);
200
    snew(ir->simtempvals, 1);
201
}
202

    
203
void stupid_fill_block(t_block *grp, int natom, gmx_bool bOneIndexGroup)
204
{
205
    int i;
206

    
207
    if (bOneIndexGroup)
208
    {
209
        grp->nalloc_index = 2;
210
        snew(grp->index, grp->nalloc_index);
211
        grp->index[0] = 0;
212
        grp->index[1] = natom;
213
        grp->nr       = 1;
214
    }
215
    else
216
    {
217
        grp->nalloc_index = natom+1;
218
        snew(grp->index, grp->nalloc_index);
219
        snew(grp->index, natom+1);
220
        for (i = 0; (i <= natom); i++)
221
        {
222
            grp->index[i] = i;
223
        }
224
        grp->nr = natom;
225
    }
226
}
227

    
228
void stupid_fill_blocka(t_blocka *grp, int natom)
229
{
230
    int i;
231

    
232
    grp->nalloc_a = natom;
233
    snew(grp->a, grp->nalloc_a);
234
    for (i = 0; (i < natom); i++)
235
    {
236
        grp->a[i] = i;
237
    }
238
    grp->nra = natom;
239

    
240
    grp->nalloc_index = natom + 1;
241
    snew(grp->index, grp->nalloc_index);
242
    for (i = 0; (i <= natom); i++)
243
    {
244
        grp->index[i] = i;
245
    }
246
    grp->nr = natom;
247
}
248

    
249
void copy_blocka(const t_blocka *src, t_blocka *dest)
250
{
251
    int i;
252

    
253
    dest->nr           = src->nr;
254
    dest->nalloc_index = dest->nr + 1;
255
    snew(dest->index, dest->nalloc_index);
256
    for (i = 0; i < dest->nr+1; i++)
257
    {
258
        dest->index[i] = src->index[i];
259
    }
260
    dest->nra      = src->nra;
261
    dest->nalloc_a = dest->nra + 1;
262
    snew(dest->a, dest->nalloc_a);
263
    for (i = 0; i < dest->nra+1; i++)
264
    {
265
        dest->a[i] = src->a[i];
266
    }
267
}
268

    
269
void done_block(t_block *block)
270
{
271
    block->nr    = 0;
272
    sfree(block->index);
273
    block->nalloc_index = 0;
274
}
275

    
276
void done_blocka(t_blocka *block)
277
{
278
    block->nr    = 0;
279
    block->nra   = 0;
280
    sfree(block->index);
281
    sfree(block->a);
282
    block->index        = NULL;
283
    block->a            = NULL;
284
    block->nalloc_index = 0;
285
    block->nalloc_a     = 0;
286
}
287

    
288
void done_atom (t_atoms *at)
289
{
290
    at->nr       = 0;
291
    at->nres     = 0;
292
    sfree(at->atom);
293
    sfree(at->resinfo);
294
    sfree(at->atomname);
295
    sfree(at->atomtype);
296
    sfree(at->atomtypeB);
297
    if (at->pdbinfo)
298
    {
299
        sfree(at->pdbinfo);
300
    }
301
}
302

    
303
void done_atomtypes(t_atomtypes *atype)
304
{
305
    atype->nr = 0;
306
    sfree(atype->radius);
307
    sfree(atype->vol);
308
    sfree(atype->surftens);
309
    sfree(atype->atomnumber);
310
    sfree(atype->gb_radius);
311
    sfree(atype->S_hct);
312
}
313

    
314
void done_moltype(gmx_moltype_t *molt)
315
{
316
    int f;
317

    
318
    done_atom(&molt->atoms);
319
    done_block(&molt->cgs);
320
    done_blocka(&molt->excls);
321

    
322
    for (f = 0; f < F_NRE; f++)
323
    {
324
        sfree(molt->ilist[f].iatoms);
325
        molt->ilist[f].nalloc = 0;
326
    }
327
}
328

    
329
void done_molblock(gmx_molblock_t *molb)
330
{
331
    if (molb->nposres_xA > 0)
332
    {
333
        molb->nposres_xA = 0;
334
        free(molb->posres_xA);
335
    }
336
    if (molb->nposres_xB > 0)
337
    {
338
        molb->nposres_xB = 0;
339
        free(molb->posres_xB);
340
    }
341
}
342

    
343
void done_mtop(gmx_mtop_t *mtop, gmx_bool bDoneSymtab)
344
{
345
    int i;
346

    
347
    if (bDoneSymtab)
348
    {
349
        done_symtab(&mtop->symtab);
350
    }
351

    
352
    sfree(mtop->ffparams.functype);
353
    sfree(mtop->ffparams.iparams);
354

    
355
    for (i = 0; i < mtop->nmoltype; i++)
356
    {
357
        done_moltype(&mtop->moltype[i]);
358
    }
359
    sfree(mtop->moltype);
360
    for (i = 0; i < mtop->nmolblock; i++)
361
    {
362
        done_molblock(&mtop->molblock[i]);
363
    }
364
    sfree(mtop->molblock);
365
    done_block(&mtop->mols);
366
}
367

    
368
void done_top(t_topology *top)
369
{
370
    int f;
371

    
372
    sfree(top->idef.functype);
373
    sfree(top->idef.iparams);
374
    for (f = 0; f < F_NRE; ++f)
375
    {
376
        sfree(top->idef.il[f].iatoms);
377
        top->idef.il[f].iatoms = NULL;
378
        top->idef.il[f].nalloc = 0;
379
    }
380

    
381
    done_atom (&(top->atoms));
382

    
383
    /* For GB */
384
    done_atomtypes(&(top->atomtypes));
385

    
386
    done_symtab(&(top->symtab));
387
    done_block(&(top->cgs));
388
    done_block(&(top->mols));
389
    done_blocka(&(top->excls));
390
}
391

    
392
static void done_pull_group(t_pull_group *pgrp)
393
{
394
    if (pgrp->nat > 0)
395
    {
396
        sfree(pgrp->ind);
397
        sfree(pgrp->ind_loc);
398
        sfree(pgrp->weight);
399
        sfree(pgrp->weight_loc);
400
    }
401
}
402

    
403
static void done_pull(t_pull *pull)
404
{
405
    int i;
406

    
407
    for (i = 0; i < pull->ngroup+1; i++)
408
    {
409
        done_pull_group(pull->group);
410
        done_pull_group(pull->dyna);
411
    }
412
}
413

    
414
void done_inputrec(t_inputrec *ir)
415
{
416
    int m;
417

    
418
    for (m = 0; (m < DIM); m++)
419
    {
420
        if (ir->ex[m].a)
421
        {
422
            sfree(ir->ex[m].a);
423
        }
424
        if (ir->ex[m].phi)
425
        {
426
            sfree(ir->ex[m].phi);
427
        }
428
        if (ir->et[m].a)
429
        {
430
            sfree(ir->et[m].a);
431
        }
432
        if (ir->et[m].phi)
433
        {
434
            sfree(ir->et[m].phi);
435
        }
436
    }
437

    
438
    sfree(ir->opts.nrdf);
439
    sfree(ir->opts.ref_t);
440
    sfree(ir->opts.annealing);
441
    sfree(ir->opts.anneal_npoints);
442
    sfree(ir->opts.anneal_time);
443
    sfree(ir->opts.anneal_temp);
444
    sfree(ir->opts.tau_t);
445
    sfree(ir->opts.acc);
446
    sfree(ir->opts.nFreeze);
447
    sfree(ir->opts.QMmethod);
448
    sfree(ir->opts.QMbasis);
449
    sfree(ir->opts.QMcharge);
450
    sfree(ir->opts.QMmult);
451
    sfree(ir->opts.bSH);
452
    sfree(ir->opts.CASorbitals);
453
    sfree(ir->opts.CASelectrons);
454
    sfree(ir->opts.SAon);
455
    sfree(ir->opts.SAoff);
456
    sfree(ir->opts.SAsteps);
457
    sfree(ir->opts.bOPT);
458
    sfree(ir->opts.bTS);
459

    
460
    if (ir->pull)
461
    {
462
        done_pull(ir->pull);
463
        sfree(ir->pull);
464
    }
465
}
466

    
467
static void zero_history(history_t *hist)
468
{
469
    hist->disre_initf  = 0;
470
    hist->ndisrepairs  = 0;
471
    hist->disre_rm3tav = NULL;
472
    hist->orire_initf  = 0;
473
    hist->norire_Dtav  = 0;
474
    hist->orire_Dtav   = NULL;
475
}
476

    
477
static void zero_ekinstate(ekinstate_t *eks)
478
{
479
    eks->ekin_n         = 0;
480
    eks->ekinh          = NULL;
481
    eks->ekinf          = NULL;
482
    eks->ekinh_old      = NULL;
483
    eks->ekinscalef_nhc = NULL;
484
    eks->ekinscaleh_nhc = NULL;
485
    eks->vscale_nhc     = NULL;
486
    eks->dekindl        = 0;
487
    eks->mvcos          = 0;
488
}
489

    
490
static void init_swapstate(swapstate_t *swapstate)
491
{
492
    int ii, ic;
493

    
494
    swapstate->eSwapCoords = 0;
495
    swapstate->nAverage    = 0;
496

    
497
    /* Ion/water position swapping */
498
    for (ic = 0; ic < eCompNR; ic++)
499
    {
500
        for (ii = 0; ii < eIonNR; ii++)
501
        {
502
            swapstate->nat_req[ic][ii]        = 0;
503
            swapstate->nat_req_p[ic][ii]      = NULL;
504
            swapstate->inflow_netto[ic][ii]   = 0;
505
            swapstate->inflow_netto_p[ic][ii] = NULL;
506
            swapstate->nat_past[ic][ii]       = NULL;
507
            swapstate->nat_past_p[ic][ii]     = NULL;
508
            swapstate->fluxfromAtoB[ic][ii]   = 0;
509
            swapstate->fluxfromAtoB_p[ic][ii] = NULL;
510
        }
511
    }
512
    swapstate->fluxleak               = NULL;
513
    swapstate->nions                  = 0;
514
    swapstate->comp_from              = NULL;
515
    swapstate->channel_label          = NULL;
516
    swapstate->bFromCpt               = 0;
517
    swapstate->nat[eChan0]            = 0;
518
    swapstate->nat[eChan1]            = 0;
519
    swapstate->xc_old_whole[eChan0]   = NULL;
520
    swapstate->xc_old_whole[eChan1]   = NULL;
521
    swapstate->xc_old_whole_p[eChan0] = NULL;
522
    swapstate->xc_old_whole_p[eChan1] = NULL;
523
}
524

    
525
void init_energyhistory(energyhistory_t * enerhist)
526
{
527
    enerhist->nener = 0;
528

    
529
    enerhist->ener_ave     = NULL;
530
    enerhist->ener_sum     = NULL;
531
    enerhist->ener_sum_sim = NULL;
532
    enerhist->dht          = NULL;
533

    
534
    enerhist->nsteps     = 0;
535
    enerhist->nsum       = 0;
536
    enerhist->nsteps_sim = 0;
537
    enerhist->nsum_sim   = 0;
538

    
539
    enerhist->dht = NULL;
540
}
541

    
542
static void done_delta_h_history(delta_h_history_t *dht)
543
{
544
    int i;
545

    
546
    for (i = 0; i < dht->nndh; i++)
547
    {
548
        sfree(dht->dh[i]);
549
    }
550
    sfree(dht->dh);
551
    sfree(dht->ndh);
552
}
553

    
554
void done_energyhistory(energyhistory_t * enerhist)
555
{
556
    sfree(enerhist->ener_ave);
557
    sfree(enerhist->ener_sum);
558
    sfree(enerhist->ener_sum_sim);
559

    
560
    if (enerhist->dht != NULL)
561
    {
562
        done_delta_h_history(enerhist->dht);
563
        sfree(enerhist->dht);
564
    }
565
}
566

    
567
void init_gtc_state(t_state *state, int ngtc, int nnhpres, int nhchainlength)
568
{
569
    int i, j;
570

    
571
    state->ngtc          = ngtc;
572
    state->nnhpres       = nnhpres;
573
    state->nhchainlength = nhchainlength;
574
    if (state->ngtc > 0)
575
    {
576
        snew(state->nosehoover_xi, state->nhchainlength*state->ngtc);
577
        snew(state->nosehoover_vxi, state->nhchainlength*state->ngtc);
578
        snew(state->therm_integral, state->ngtc);
579
        for (i = 0; i < state->ngtc; i++)
580
        {
581
            for (j = 0; j < state->nhchainlength; j++)
582
            {
583
                state->nosehoover_xi[i*state->nhchainlength + j]   = 0.0;
584
                state->nosehoover_vxi[i*state->nhchainlength + j]  = 0.0;
585
            }
586
        }
587
        for (i = 0; i < state->ngtc; i++)
588
        {
589
            state->therm_integral[i]  = 0.0;
590
        }
591
    }
592
    else
593
    {
594
        state->nosehoover_xi  = NULL;
595
        state->nosehoover_vxi = NULL;
596
        state->therm_integral = NULL;
597
    }
598

    
599
    if (state->nnhpres > 0)
600
    {
601
        snew(state->nhpres_xi, state->nhchainlength*nnhpres);
602
        snew(state->nhpres_vxi, state->nhchainlength*nnhpres);
603
        for (i = 0; i < nnhpres; i++)
604
        {
605
            for (j = 0; j < state->nhchainlength; j++)
606
            {
607
                state->nhpres_xi[i*nhchainlength + j]   = 0.0;
608
                state->nhpres_vxi[i*nhchainlength + j]  = 0.0;
609
            }
610
        }
611
    }
612
    else
613
    {
614
        state->nhpres_xi  = NULL;
615
        state->nhpres_vxi = NULL;
616
    }
617
}
618

    
619

    
620
void init_state(t_state *state, int natoms, int ngtc, int nnhpres, int nhchainlength, int nlambda)
621
{
622
    int i;
623

    
624
    state->natoms = natoms;
625
    state->flags  = 0;
626
    state->lambda = 0;
627
    snew(state->lambda, efptNR);
628
    for (i = 0; i < efptNR; i++)
629
    {
630
        state->lambda[i] = 0;
631
    }
632
    state->veta   = 0;
633
    clear_mat(state->box);
634
    clear_mat(state->box_rel);
635
    clear_mat(state->boxv);
636
    clear_mat(state->pres_prev);
637
    clear_mat(state->svir_prev);
638
    clear_mat(state->fvir_prev);
639
    init_gtc_state(state, ngtc, nnhpres, nhchainlength);
640
    state->nalloc = state->natoms;
641
    if (state->nalloc > 0)
642
    {
643
        snew(state->x, state->nalloc);
644
        snew(state->v, state->nalloc);
645
    }
646
    else
647
    {
648
        state->x = NULL;
649
        state->v = NULL;
650
    }
651
    state->sd_X = NULL;
652
    state->cg_p = NULL;
653
    zero_history(&state->hist);
654
    zero_ekinstate(&state->ekinstate);
655
    init_energyhistory(&state->enerhist);
656
    init_df_history(&state->dfhist, nlambda);
657
    init_swapstate(&state->swapstate);
658
    state->ddp_count       = 0;
659
    state->ddp_count_cg_gl = 0;
660
    state->cg_gl           = NULL;
661
    state->cg_gl_nalloc    = 0;
662
}
663

    
664
void done_state(t_state *state)
665
{
666
    if (state->x)
667
    {
668
        sfree(state->x);
669
    }
670
    if (state->v)
671
    {
672
        sfree(state->v);
673
    }
674
    if (state->sd_X)
675
    {
676
        sfree(state->sd_X);
677
    }
678
    if (state->cg_p)
679
    {
680
        sfree(state->cg_p);
681
    }
682
    state->nalloc = 0;
683
    if (state->cg_gl)
684
    {
685
        sfree(state->cg_gl);
686
    }
687
    state->cg_gl_nalloc = 0;
688
    if (state->lambda)
689
    {
690
        sfree(state->lambda);
691
    }
692
    if (state->ngtc > 0)
693
    {
694
        sfree(state->nosehoover_xi);
695
        sfree(state->nosehoover_vxi);
696
        sfree(state->therm_integral);
697
    }
698
}
699

    
700
t_state *serial_init_local_state(t_state *state_global)
701
{
702
    int      i;
703
    t_state *state_local;
704

    
705
    snew(state_local, 1);
706

    
707
    /* Copy all the contents */
708
    *state_local = *state_global;
709
    snew(state_local->lambda, efptNR);
710
    /* local storage for lambda */
711
    for (i = 0; i < efptNR; i++)
712
    {
713
        state_local->lambda[i] = state_global->lambda[i];
714
    }
715

    
716
    return state_local;
717
}
718

    
719
static void do_box_rel(t_inputrec *ir, matrix box_rel, matrix b, gmx_bool bInit)
720
{
721
    int d, d2;
722

    
723
    for (d = YY; d <= ZZ; d++)
724
    {
725
        for (d2 = XX; d2 <= (ir->epct == epctSEMIISOTROPIC ? YY : ZZ); d2++)
726
        {
727
            /* We need to check if this box component is deformed
728
             * or if deformation of another component might cause
729
             * changes in this component due to box corrections.
730
             */
731
            if (ir->deform[d][d2] == 0 &&
732
                !(d == ZZ && d2 == XX && ir->deform[d][YY] != 0 &&
733
                  (b[YY][d2] != 0 || ir->deform[YY][d2] != 0)))
734
            {
735
                if (bInit)
736
                {
737
                    box_rel[d][d2] = b[d][d2]/b[XX][XX];
738
                }
739
                else
740
                {
741
                    b[d][d2] = b[XX][XX]*box_rel[d][d2];
742
                }
743
            }
744
        }
745
    }
746
}
747

    
748
void set_box_rel(t_inputrec *ir, t_state *state)
749
{
750
    /* Make sure the box obeys the restrictions before we fix the ratios */
751
    correct_box(NULL, 0, state->box, NULL);
752

    
753
    clear_mat(state->box_rel);
754

    
755
    if (PRESERVE_SHAPE(*ir))
756
    {
757
        do_box_rel(ir, state->box_rel, state->box, TRUE);
758
    }
759
}
760

    
761
void preserve_box_shape(t_inputrec *ir, matrix box_rel, matrix b)
762
{
763
    if (PRESERVE_SHAPE(*ir))
764
    {
765
        do_box_rel(ir, box_rel, b, FALSE);
766
    }
767
}
768

    
769
void add_t_atoms(t_atoms *atoms, int natom_extra, int nres_extra)
770
{
771
    int i;
772

    
773
    if (natom_extra > 0)
774
    {
775
        srenew(atoms->atomname, atoms->nr+natom_extra);
776
        srenew(atoms->atom, atoms->nr+natom_extra);
777
        if (NULL != atoms->pdbinfo)
778
        {
779
            srenew(atoms->pdbinfo, atoms->nr+natom_extra);
780
        }
781
        if (NULL != atoms->atomtype)
782
        {
783
            srenew(atoms->atomtype, atoms->nr+natom_extra);
784
        }
785
        if (NULL != atoms->atomtypeB)
786
        {
787
            srenew(atoms->atomtypeB, atoms->nr+natom_extra);
788
        }
789
        for (i = atoms->nr; (i < atoms->nr+natom_extra); i++)
790
        {
791
            atoms->atomname[i] = NULL;
792
            memset(&atoms->atom[i], 0, sizeof(atoms->atom[i]));
793
            if (NULL != atoms->pdbinfo)
794
            {
795
                memset(&atoms->pdbinfo[i], 0, sizeof(atoms->pdbinfo[i]));
796
            }
797
            if (NULL != atoms->atomtype)
798
            {
799
                atoms->atomtype[i] = NULL;
800
            }
801
            if (NULL != atoms->atomtypeB)
802
            {
803
                atoms->atomtypeB[i] = NULL;
804
            }
805
        }
806
        atoms->nr += natom_extra;
807
    }
808
    if (nres_extra > 0)
809
    {
810
        srenew(atoms->resinfo, atoms->nres+nres_extra);
811
        for (i = atoms->nres; (i < atoms->nres+nres_extra); i++)
812
        {
813
            memset(&atoms->resinfo[i], 0, sizeof(atoms->resinfo[i]));
814
        }
815
        atoms->nres += nres_extra;
816
    }
817
}
818

    
819
void init_t_atoms(t_atoms *atoms, int natoms, gmx_bool bPdbinfo)
820
{
821
    atoms->nr   = natoms;
822
    atoms->nres = 0;
823
    snew(atoms->atomname, natoms);
824
    atoms->atomtype  = NULL;
825
    atoms->atomtypeB = NULL;
826
    snew(atoms->resinfo, natoms);
827
    snew(atoms->atom, natoms);
828
    if (bPdbinfo)
829
    {
830
        snew(atoms->pdbinfo, natoms);
831
    }
832
    else
833
    {
834
        atoms->pdbinfo = NULL;
835
    }
836
}
837

    
838
t_atoms *copy_t_atoms(t_atoms *src)
839
{
840
    t_atoms *dst;
841
    int      i;
842

    
843
    snew(dst, 1);
844
    init_t_atoms(dst, src->nr, (NULL != src->pdbinfo));
845
    dst->nr = src->nr;
846
    if (NULL != src->atomname)
847
    {
848
        snew(dst->atomname, src->nr);
849
    }
850
    if (NULL != src->atomtype)
851
    {
852
        snew(dst->atomtype, src->nr);
853
    }
854
    if (NULL != src->atomtypeB)
855
    {
856
        snew(dst->atomtypeB, src->nr);
857
    }
858
    for (i = 0; (i < src->nr); i++)
859
    {
860
        dst->atom[i] = src->atom[i];
861
        if (NULL != src->pdbinfo)
862
        {
863
            dst->pdbinfo[i] = src->pdbinfo[i];
864
        }
865
        if (NULL != src->atomname)
866
        {
867
            dst->atomname[i]  = src->atomname[i];
868
        }
869
        if (NULL != src->atomtype)
870
        {
871
            dst->atomtype[i] = src->atomtype[i];
872
        }
873
        if (NULL != src->atomtypeB)
874
        {
875
            dst->atomtypeB[i] = src->atomtypeB[i];
876
        }
877
    }
878
    dst->nres = src->nres;
879
    for (i = 0; (i < src->nres); i++)
880
    {
881
        dst->resinfo[i] = src->resinfo[i];
882
    }
883
    return dst;
884
}
885

    
886
void t_atoms_set_resinfo(t_atoms *atoms, int atom_ind, t_symtab *symtab,
887
                         const char *resname, int resnr, unsigned char ic,
888
                         int chainnum, char chainid)
889
{
890
    t_resinfo *ri;
891

    
892
    ri           = &atoms->resinfo[atoms->atom[atom_ind].resind];
893
    ri->name     = put_symtab(symtab, resname);
894
    ri->rtp      = NULL;
895
    ri->nr       = resnr;
896
    ri->ic       = ic;
897
    ri->chainnum = chainnum;
898
    ri->chainid  = chainid;
899
}
900

    
901
void free_t_atoms(t_atoms *atoms, gmx_bool bFreeNames)
902
{
903
    int i;
904

    
905
    if (bFreeNames && atoms->atomname != NULL)
906
    {
907
        for (i = 0; i < atoms->nr; i++)
908
        {
909
            if (atoms->atomname[i] != NULL)
910
            {
911
                sfree(*atoms->atomname[i]);
912
                *atoms->atomname[i] = NULL;
913
            }
914
        }
915
    }
916
    if (bFreeNames && atoms->resinfo != NULL)
917
    {
918
        for (i = 0; i < atoms->nres; i++)
919
        {
920
            if (atoms->resinfo[i].name != NULL)
921
            {
922
                sfree(*atoms->resinfo[i].name);
923
                *atoms->resinfo[i].name = NULL;
924
            }
925
        }
926
    }
927
    if (bFreeNames && atoms->atomtype != NULL)
928
    {
929
        for (i = 0; i < atoms->nr; i++)
930
        {
931
            if (atoms->atomtype[i] != NULL)
932
            {
933
                sfree(*atoms->atomtype[i]);
934
                *atoms->atomtype[i] = NULL;
935
            }
936
        }
937
    }
938
    if (bFreeNames && atoms->atomtypeB != NULL)
939
    {
940
        for (i = 0; i < atoms->nr; i++)
941
        {
942
            if (atoms->atomtypeB[i] != NULL)
943
            {
944
                sfree(*atoms->atomtypeB[i]);
945
                *atoms->atomtypeB[i] = NULL;
946
            }
947
        }
948
    }
949
    sfree(atoms->atomname);
950
    sfree(atoms->atomtype);
951
    sfree(atoms->atomtypeB);
952
    sfree(atoms->resinfo);
953
    sfree(atoms->atom);
954
    sfree(atoms->pdbinfo);
955
    atoms->nr        = 0;
956
    atoms->nres      = 0;
957
    atoms->atomname  = NULL;
958
    atoms->atomtype  = NULL;
959
    atoms->atomtypeB = NULL;
960
    atoms->resinfo   = NULL;
961
    atoms->atom      = NULL;
962
    atoms->pdbinfo   = NULL;
963
}
964

    
965
real max_cutoff(real cutoff1, real cutoff2)
966
{
967
    if (cutoff1 == 0 || cutoff2 == 0)
968
    {
969
        return 0;
970
    }
971
    else
972
    {
973
        return max(cutoff1, cutoff2);
974
    }
975
}
976

    
977
void init_df_history(df_history_t *dfhist, int nlambda)
978
{
979
    int i;
980

    
981
    dfhist->nlambda  = nlambda;
982
    dfhist->bEquil   = 0;
983
    dfhist->wl_delta = 0;
984

    
985
    if (nlambda > 0)
986
    {
987
        snew(dfhist->sum_weights, dfhist->nlambda);
988
        snew(dfhist->sum_dg, dfhist->nlambda);
989
        snew(dfhist->sum_minvar, dfhist->nlambda);
990
        snew(dfhist->sum_variance, dfhist->nlambda);
991
        snew(dfhist->n_at_lam, dfhist->nlambda);
992
        snew(dfhist->wl_histo, dfhist->nlambda);
993
        snew(dfhist->dfavg, dfhist->nlambda);
994
        snew(dfhist->sum_dvdl, dfhist->nlambda);
995
        snew(dfhist->store_fepot, dfhist->nlambda);
996
        snew(dfhist->aim_at_lam, dfhist->nlambda);
997
    snew(dfhist->laccept, dfhist->nlambda);
998

    
999
        /* allocate transition matrices here */
1000
        snew(dfhist->Tij, dfhist->nlambda);
1001
        snew(dfhist->Tij_empirical, dfhist->nlambda);
1002

    
1003
        /* allocate accumulators for various transition matrix
1004
           free energy methods here */
1005
        snew(dfhist->accum_p, dfhist->nlambda);
1006
        snew(dfhist->accum_m, dfhist->nlambda);
1007
        snew(dfhist->accum_p2, dfhist->nlambda);
1008
        snew(dfhist->accum_m2, dfhist->nlambda);
1009

    
1010
        for (i = 0; i < dfhist->nlambda; i++)
1011
        {
1012
            snew(dfhist->Tij[i], dfhist->nlambda);
1013
            snew(dfhist->Tij_empirical[i], dfhist->nlambda);
1014
            snew((dfhist->accum_p)[i], dfhist->nlambda);
1015
            snew((dfhist->accum_m)[i], dfhist->nlambda);
1016
            snew((dfhist->accum_p2)[i], dfhist->nlambda);
1017
            snew((dfhist->accum_m2)[i], dfhist->nlambda);
1018
        }
1019
    }
1020
}
1021

    
1022
extern void copy_df_history(df_history_t *df_dest, df_history_t *df_source)
1023
{
1024
    int i, j;
1025

    
1026
    /* Currently, there should not be any difference in nlambda between the two,
1027
       but this is included for completeness for potential later functionality */
1028
    df_dest->nlambda  = df_source->nlambda;
1029
    df_dest->bEquil   = df_source->bEquil;
1030
    df_dest->wl_delta = df_source->wl_delta;
1031

    
1032
    for (i = 0; i < df_dest->nlambda; i++)
1033
    {
1034
        df_dest->sum_weights[i]  = df_source->sum_weights[i];
1035
        df_dest->sum_dg[i]       = df_source->sum_dg[i];
1036
        df_dest->sum_minvar[i]   = df_source->sum_minvar[i];
1037
        df_dest->sum_variance[i] = df_source->sum_variance[i];
1038
        df_dest->n_at_lam[i]     = df_source->n_at_lam[i];
1039
        df_dest->wl_histo[i]     = df_source->wl_histo[i];
1040
        df_dest->dfavg[i]        = df_source->dfavg[i];
1041
        df_dest->sum_dvdl[i]     = df_source->sum_dvdl[i];
1042
        df_dest->store_fepot[i]  = df_source->store_fepot[i];
1043
        df_dest->aim_at_lam[i]   = df_source->aim_at_lam[i];
1044
    df_dest->laccept[i]      = df_source->laccept[i];
1045
    }
1046

    
1047
    for (i = 0; i < df_dest->nlambda; i++)
1048
    {
1049
        for (j = 0; j < df_dest->nlambda; j++)
1050
        {
1051
            df_dest->accum_p[i][j]        = df_source->accum_p[i][j];
1052
            df_dest->accum_m[i][j]        = df_source->accum_m[i][j];
1053
            df_dest->accum_p2[i][j]       = df_source->accum_p2[i][j];
1054
            df_dest->accum_m2[i][j]       = df_source->accum_m2[i][j];
1055
            df_dest->Tij[i][j]            = df_source->Tij[i][j];
1056
            df_dest->Tij_empirical[i][j]  = df_source->Tij_empirical[i][j];
1057
        }
1058
    }
1059
}
1060

    
1061
void done_df_history(df_history_t *dfhist)
1062
{
1063
    int i;
1064

    
1065
    if (dfhist->nlambda > 0)
1066
    {
1067
        sfree(dfhist->n_at_lam);
1068
        sfree(dfhist->wl_histo);
1069
        sfree(dfhist->sum_weights);
1070
        sfree(dfhist->sum_dg);
1071
        sfree(dfhist->sum_minvar);
1072
        sfree(dfhist->sum_variance);
1073
        sfree(dfhist->dfavg);
1074
        sfree(dfhist->sum_dvdl);
1075
        sfree(dfhist->store_fepot);
1076
        sfree(dfhist->aim_at_lam);
1077
    sfree(dfhist->laccept);
1078

    
1079
        for (i = 0; i < dfhist->nlambda; i++)
1080
        {
1081
            sfree(dfhist->Tij[i]);
1082
            sfree(dfhist->Tij_empirical[i]);
1083
            sfree(dfhist->accum_p[i]);
1084
            sfree(dfhist->accum_m[i]);
1085
            sfree(dfhist->accum_p2[i]);
1086
            sfree(dfhist->accum_m2[i]);
1087
        }
1088
    }
1089
    dfhist->bEquil   = 0;
1090
    dfhist->nlambda  = 0;
1091
    dfhist->wl_delta = 0;
1092
}