Project

General

Profile

gmx_density.c

klark chen, 09/05/2017 09:43 AM

 
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,2015, 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
#include "gmxpre.h"
38

    
39
#include <ctype.h>
40
#include <math.h>
41
#include <stdlib.h>
42
#include <string.h>
43

    
44
#include "gromacs/commandline/pargs.h"
45
#include "gromacs/fileio/tpxio.h"
46
#include "gromacs/fileio/trxio.h"
47
#include "gromacs/fileio/xvgr.h"
48
#include "gromacs/gmxana/gmx_ana.h"
49
#include "gromacs/gmxana/gstat.h"
50
#include "gromacs/legacyheaders/macros.h"
51
#include "gromacs/legacyheaders/typedefs.h"
52
#include "gromacs/legacyheaders/viewit.h"
53
#include "gromacs/math/units.h"
54
#include "gromacs/math/vec.h"
55
#include "gromacs/pbcutil/pbc.h"
56
#include "gromacs/pbcutil/rmpbc.h"
57
#include "gromacs/topology/index.h"
58
#include "gromacs/utility/cstringutil.h"
59
#include "gromacs/utility/fatalerror.h"
60
#include "gromacs/utility/futil.h"
61
#include "gromacs/utility/smalloc.h"
62

    
63
typedef struct {
64
    char *atomname;
65
    int   nr_el;
66
} t_electron;
67

    
68
/****************************************************************************/
69
/* This program calculates the partial density across the box.              */
70
/* Peter Tieleman, Mei 1995                                                 */
71
/****************************************************************************/
72

    
73
/* used for sorting the list */
74
int compare(void *a, void *b)
75
{
76
    t_electron *tmp1, *tmp2;
77
    tmp1 = (t_electron *)a; tmp2 = (t_electron *)b;
78

    
79
    return strcmp(tmp1->atomname, tmp2->atomname);
80
}
81

    
82
int get_electrons(t_electron **eltab, const char *fn)
83
{
84
    char  buffer[256];  /* to read in a line   */
85
    char  tempname[80]; /* buffer to hold name */
86
    int   tempnr;
87

    
88
    FILE *in;
89
    int   nr;        /* number of atomstypes to read */
90
    int   i;
91

    
92
    if (!(in = gmx_ffopen(fn, "r")))
93
    {
94
        gmx_fatal(FARGS, "Couldn't open %s. Exiting.\n", fn);
95
    }
96

    
97
    if (NULL == fgets(buffer, 255, in))
98
    {
99
        gmx_fatal(FARGS, "Error reading from file %s", fn);
100
    }
101

    
102
    if (sscanf(buffer, "%d", &nr) != 1)
103
    {
104
        gmx_fatal(FARGS, "Invalid number of atomtypes in datafile\n");
105
    }
106

    
107
    snew(*eltab, nr);
108

    
109
    for (i = 0; i < nr; i++)
110
    {
111
        if (fgets(buffer, 255, in) == NULL)
112
        {
113
            gmx_fatal(FARGS, "reading datafile. Check your datafile.\n");
114
        }
115
        if (sscanf(buffer, "%s = %d", tempname, &tempnr) != 2)
116
        {
117
            gmx_fatal(FARGS, "Invalid line in datafile at line %d\n", i+1);
118
        }
119
        (*eltab)[i].nr_el    = tempnr;
120
        (*eltab)[i].atomname = gmx_strdup(tempname);
121
    }
122
    gmx_ffclose(in);
123

    
124
    /* sort the list */
125
    fprintf(stderr, "Sorting list..\n");
126
    qsort ((void*)*eltab, nr, sizeof(t_electron),
127
           (int(*)(const void*, const void*))compare);
128

    
129
    return nr;
130
}
131

    
132
void center_coords(t_atoms *atoms, atom_id *index_center, int ncenter,
133
                   matrix box, rvec x0[])
134
{
135
    int  i, k, m;
136
    real tmass, mm;
137
    rvec com, shift, box_center;
138

    
139
    tmass = 0;
140
    clear_rvec(com);
141
    for (k = 0; (k < ncenter); k++)
142
    {
143
        i = index_center[k];
144
        if (i >= atoms->nr)
145
        {
146
            gmx_fatal(FARGS, "Index %d refers to atom %d, which is larger than natoms (%d).",
147
                      k+1, i+1, atoms->nr);
148
        }
149
        mm     = atoms->atom[i].m;
150
        tmass += mm;
151
        for (m = 0; (m < DIM); m++)
152
        {
153
            com[m] += mm*x0[i][m];
154
        }
155
    }
156
    for (m = 0; (m < DIM); m++)
157
    {
158
        com[m] /= tmass;
159
    }
160
    calc_box_center(ecenterDEF, box, box_center);
161
    rvec_sub(com, box_center, shift);
162

    
163
    /* Important - while the center was calculated based on a group, we should move all atoms */
164
    for (i = 0; (i < atoms->nr); i++)
165
    {
166
        rvec_dec(x0[i], shift);
167
    }
168
}
169

    
170
void calc_electron_density(const char *fn, atom_id **index, int gnx[],
171
                           double ***slDensity, int *nslices, t_topology *top,
172
                           int ePBC,
173
                           int axis, int nr_grps, real *slWidth,
174
                           t_electron eltab[], int nr, gmx_bool bCenter,
175
                           atom_id *index_center, int ncenter,
176
                           gmx_bool bRelative, const output_env_t oenv)
177
{
178
    rvec        *x0;            /* coordinates without pbc */
179
    matrix       box;           /* box (3x3) */
180
    double       invvol;
181
    int          natoms;        /* nr. atoms in trj */
182
    t_trxstatus *status;
183
    int          i, n,          /* loop indices */
184
                 nr_frames = 0, /* number of frames */
185
                 slice;         /* current slice */
186
    t_electron  *found;         /* found by bsearch */
187
    t_electron   sought;        /* thingie thought by bsearch */
188
    real         boxSz, aveBox;
189
    gmx_rmpbc_t  gpbc = NULL;
190

    
191
    real         t,
192
                 z;
193

    
194
    if (axis < 0 || axis >= DIM)
195
    {
196
        gmx_fatal(FARGS, "Invalid axes. Terminating\n");
197
    }
198

    
199
    if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
200
    {
201
        gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
202
    }
203

    
204
    aveBox = 0;
205

    
206
    if (!*nslices)
207
    {
208
        *nslices = (int)(box[axis][axis] * 10); /* default value */
209
        fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
210
    }
211

    
212
    snew(*slDensity, nr_grps);
213
    for (i = 0; i < nr_grps; i++)
214
    {
215
        snew((*slDensity)[i], *nslices);
216
    }
217

    
218
    gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
219
    /*********** Start processing trajectory ***********/
220
    do
221
    {
222
        gmx_rmpbc(gpbc, natoms, box, x0);
223

    
224
        /* Translate atoms so the com of the center-group is in the
225
         * box geometrical center.
226
         */
227
        if (bCenter)
228
        {
229
            center_coords(&top->atoms, index_center, ncenter, box, x0);
230
        }
231

    
232
        invvol   = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
233

    
234
        if (bRelative)
235
        {
236
            *slWidth = 1.0/(*nslices);
237
            boxSz    = 1.0;
238
        }
239
        else
240
        {
241
            *slWidth = box[axis][axis]/(*nslices);
242
            boxSz    = box[axis][axis];
243
        }
244

    
245
        aveBox += box[axis][axis];
246

    
247
        for (n = 0; n < nr_grps; n++)
248
        {
249
            for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
250
            {
251
                z = x0[index[n][i]][axis];
252
                while (z < 0)
253
                {
254
                    z += box[axis][axis];
255
                }
256
                while (z > box[axis][axis])
257
                {
258
                    z -= box[axis][axis];
259
                }
260

    
261
                if (bRelative)
262
                {
263
                    z = z/box[axis][axis];
264
                }
265

    
266
                /* determine which slice atom is in */
267
                if (bCenter)
268
                {
269
                    slice = floor( (z-(boxSz/2.0)) / (*slWidth) ) + *nslices/2;
270
                }
271
                else
272
                {
273
                    slice = (z / (*slWidth));
274
                }
275
                sought.nr_el    = 0;
276
                sought.atomname = gmx_strdup(*(top->atoms.atomname[index[n][i]]));
277

    
278
                /* now find the number of electrons. This is not efficient. */
279
                found = (t_electron *)
280
                    bsearch((const void *)&sought,
281
                            (const void *)eltab, nr, sizeof(t_electron),
282
                            (int(*)(const void*, const void*))compare);
283

    
284
                if (found == NULL)
285
                {
286
                    fprintf(stderr, "Couldn't find %s. Add it to the .dat file\n",
287
                            *(top->atoms.atomname[index[n][i]]));
288
                }
289
                else
290
                {
291
                    (*slDensity)[n][slice] += (found->nr_el -
292
                                               top->atoms.atom[index[n][i]].q)*invvol;
293
                }
294
                free(sought.atomname);
295
            }
296
        }
297
        nr_frames++;
298
    }
299
    while (read_next_x(oenv, status, &t, x0, box));
300
    gmx_rmpbc_done(gpbc);
301

    
302
    /*********** done with status file **********/
303
    close_trj(status);
304

    
305
/* slDensity now contains the total number of electrons per slice, summed
306
   over all frames. Now divide by nr_frames and volume of slice
307
 */
308

    
309
    fprintf(stderr, "\nRead %d frames from trajectory. Counting electrons\n",
310
            nr_frames);
311

    
312
    if (bRelative)
313
    {
314
        aveBox  /= nr_frames;
315
        *slWidth = aveBox/(*nslices);
316
    }
317

    
318
    for (n = 0; n < nr_grps; n++)
319
    {
320
        for (i = 0; i < *nslices; i++)
321
        {
322
            (*slDensity)[n][i] /= nr_frames;
323
        }
324
    }
325

    
326
    sfree(x0); /* free memory used by coordinate array */
327
}
328

    
329
void calc_density(const char *fn, atom_id **index, int gnx[],
330
                  double ***slDensity, int *nslices, t_topology *top, int ePBC,
331
                  int axis, int nr_grps, real *slWidth, gmx_bool bCenter,
332
                  atom_id *index_center, int ncenter,
333
                  gmx_bool bRelative, const output_env_t oenv, const char **dens_opt)
334
{
335
    rvec        *x0;            /* coordinates without pbc */
336
    matrix       box;           /* box (3x3) */
337
    double       invvol;
338
    int          natoms;        /* nr. atoms in trj */
339
    t_trxstatus *status;
340
    int        **slCount,       /* nr. of atoms in one slice for a group */
341
                 i, j, n,       /* loop indices */
342
                 ax1       = 0, ax2 = 0,
343
                 nr_frames = 0, /* number of frames */
344
                 slice;         /* current slice */
345
    real         t,
346
                 z;
347
    real         boxSz, aveBox;
348
    char        *buf;    /* for tmp. keeping atomname */
349
    real        *den_val; /* values from which the density is calculated */
350
    gmx_rmpbc_t  gpbc = NULL;
351

    
352
    if (axis < 0 || axis >= DIM)
353
    {
354
        gmx_fatal(FARGS, "Invalid axes. Terminating\n");
355
    }
356

    
357
    if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
358
    {
359
        gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
360
    }
361

    
362
    aveBox = 0;
363

    
364
    if (!*nslices)
365
    {
366
        *nslices = (int)(box[axis][axis] * 10); /* default value */
367
        fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
368
    }
369

    
370
    snew(*slDensity, nr_grps);
371
    for (i = 0; i < nr_grps; i++)
372
    {
373
        snew((*slDensity)[i], *nslices);
374
    }
375

    
376
    gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
377
    /*********** Start processing trajectory ***********/
378

    
379
    snew(den_val, top->atoms.nr);
380
    if (dens_opt[0][0] == 'n')
381
    {
382
        for (i = 0; (i < top->atoms.nr); i++)
383
        {
384
            den_val[i] = 1;
385
        }
386
    }
387
    else if (dens_opt[0][0] == 'c')
388
    {
389
        for (i = 0; (i < top->atoms.nr); i++)
390
        {
391
            den_val[i] = top->atoms.atom[i].q;
392
        }
393
    }
394
    else
395
    {
396
        for (i = 0; (i < top->atoms.nr); i++)
397
        {
398
            den_val[i] = top->atoms.atom[i].m;
399
        }
400
    }
401

    
402
    do
403
    {
404
        gmx_rmpbc(gpbc, natoms, box, x0);
405

    
406
        /* Translate atoms so the com of the center-group is in the
407
         * box geometrical center.
408
         */
409
        if (bCenter)
410
        {
411
            center_coords(&top->atoms, index_center, ncenter, box, x0);
412
        }
413

    
414
        invvol   = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
415

    
416
        if (bRelative)
417
        {
418
            *slWidth = 1.0/(*nslices);
419
            boxSz    = 1.0;
420
        }
421
        else
422
        {
423
            *slWidth = box[axis][axis]/(*nslices);
424
            boxSz    = box[axis][axis];
425
        }
426

    
427
        aveBox += box[axis][axis];
428

    
429
        for (n = 0; n < nr_grps; n++)
430
        {
431
            for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
432
            {
433
                z = x0[index[n][i]][axis];
434
                while (z < 0)
435
                {
436
                    z += box[axis][axis];
437
                }
438
                while (z > box[axis][axis])
439
                {
440
                    z -= box[axis][axis];
441
                }
442

    
443
                if (bRelative)
444
                {
445
                    z = z/box[axis][axis];
446
                }
447

    
448
                /* determine which slice atom is in */
449
                if (bCenter)
450
                {
451
                    slice = floor( (z-(boxSz/2.0)) / (*slWidth) ) + *nslices/2;
452
                }
453
                else
454
                {
455
                    slice = floor(z / (*slWidth));
456
                }
457

    
458
                /* Slice should already be 0<=slice<nslices, but we just make
459
                 * sure we are not hit by IEEE rounding errors since we do
460
                 * math operations after applying PBC above.
461
                 */
462
                if (slice < 0)
463
                {
464
                    slice += *nslices;
465
                }
466
                else if (slice >= *nslices)
467
                {
468
                    slice -= *nslices;
469
                }
470

    
471
                (*slDensity)[n][slice] += den_val[index[n][i]]*invvol;
472
            }
473
        }
474
        nr_frames++;
475
    }
476
    while (read_next_x(oenv, status, &t, x0, box));
477
    gmx_rmpbc_done(gpbc);
478

    
479
    /*********** done with status file **********/
480
    close_trj(status);
481

    
482
    /* slDensity now contains the total mass per slice, summed over all
483
       frames. Now divide by nr_frames and volume of slice
484
     */
485

    
486
    fprintf(stderr, "\nRead %d frames from trajectory. Calculating density\n",
487
            nr_frames);
488

    
489
    if (bRelative)
490
    {
491
        aveBox  /= nr_frames;
492
        *slWidth = aveBox/(*nslices);
493
    }
494

    
495
    for (n = 0; n < nr_grps; n++)
496
    {
497
        for (i = 0; i < *nslices; i++)
498
        {
499
            (*slDensity)[n][i] /= nr_frames;
500
        }
501
    }
502

    
503
    sfree(x0); /* free memory used by coordinate array */
504
    sfree(den_val);
505
}
506

    
507
void plot_density(double *slDensity[], const char *afile, int nslices,
508
                  int nr_grps, char *grpname[], real slWidth,
509
                  const char **dens_opt,
510
                  gmx_bool bCenter, gmx_bool bRelative, gmx_bool bSymmetrize,
511
                  const output_env_t oenv)
512
{
513
    FILE       *den;
514
    const char *title  = NULL;
515
    const char *xlabel = NULL;
516
    const char *ylabel = NULL;
517
    int         slice, n;
518
    real        ddd;
519
    real        axispos;
520

    
521
    title = bSymmetrize ? "Symmetrized partial density" : "Partial density";
522

    
523
    if (bCenter)
524
    {
525
        xlabel = bRelative ?
526
            "Average relative position from center (nm)" :
527
            "Relative position from center (nm)";
528
    }
529
    else
530
    {
531
        xlabel = bRelative ? "Average coordinate (nm)" : "Coordinate (nm)";
532
    }
533

    
534
    switch (dens_opt[0][0])
535
    {
536
        case 'm': ylabel = "Density (kg m\\S-3\\N)"; break;
537
        case 'n': ylabel = "Number density (nm\\S-3\\N)"; break;
538
        case 'c': ylabel = "Charge density (e nm\\S-3\\N)"; break;
539
        case 'e': ylabel = "Electron density (e nm\\S-3\\N)"; break;
540
    }
541

    
542
    den = xvgropen(afile,
543
                   title, xlabel, ylabel, oenv);
544

    
545
    xvgr_legend(den, nr_grps, (const char**)grpname, oenv);
546

    
547
    for (slice = 0; (slice < nslices); slice++)
548
    {
549
        if (bCenter)
550
        {
551
            axispos = (slice - nslices/2.0 + 0.5) * slWidth;
552
        }
553
        else
554
        {
555
            axispos = (slice + 0.5) * slWidth;
556
        }
557
        fprintf(den, "%12g  ", axispos);
558
        for (n = 0; (n < nr_grps); n++)
559
        {
560
            if (bSymmetrize)
561
            {
562
                ddd = (slDensity[n][slice]+slDensity[n][nslices-slice-1])*0.5;
563
            }
564
            else
565
            {
566
                ddd = slDensity[n][slice];
567
            }
568
            if (dens_opt[0][0] == 'm')
569
            {
570
                fprintf(den, "   %12g", ddd*AMU/(NANO*NANO*NANO));
571
            }
572
            else
573
            {
574
                fprintf(den, "   %12g", ddd);
575
            }
576
        }
577
        fprintf(den, "\n");
578
    }
579

    
580
    xvgrclose(den);
581
}
582

    
583
int gmx_density(int argc, char *argv[])
584
{
585
    const char        *desc[] = {
586
        "[THISMODULE] computes partial densities across the box, using an index file.[PAR]",
587
        "For the total density of NPT simulations, use [gmx-energy] instead.",
588
        "[PAR]",
589

    
590
        "Option [TT]-center[tt] performs the histogram binning relative to the center",
591
        "of an arbitrary group, in absolute box coordinates. If you are calculating",
592
        "profiles along the Z axis box dimension bZ, output would be from -bZ/2 to",
593
        "bZ/2 if you center based on the entire system.",
594
        "Note that this behaviour has changed in GROMACS 5.0; earlier versions",
595
        "merely performed a static binning in (0,bZ) and shifted the output. Now",
596
        "we compute the center for each frame and bin in (-bZ/2,bZ/2).[PAR]",
597

    
598
        "Option [TT]-symm[tt] symmetrizes the output around the center. This will",
599
        "automatically turn on [TT]-center[tt] too.",
600

    
601
        "Option [TT]-relative[tt] performs the binning in relative instead of absolute",
602
        "box coordinates, and scales the final output with the average box dimension",
603
        "along the output axis. This can be used in combination with [TT]-center[tt].[PAR]",
604

    
605
        "Densities are in kg/m^3, and number densities or electron densities can also be",
606
        "calculated. For electron densities, a file describing the number of",
607
        "electrons for each type of atom should be provided using [TT]-ei[tt].",
608
        "It should look like::",
609
        "",
610
        "   2",
611
        "   atomname = nrelectrons",
612
        "   atomname = nrelectrons",
613
        "",
614
        "The first line contains the number of lines to read from the file.",
615
        "There should be one line for each unique atom name in your system.",
616
        "The number of electrons for each atom is modified by its atomic",
617
        "partial charge.[PAR]",
618

    
619
        "IMPORTANT CONSIDERATIONS FOR BILAYERS[PAR]",
620
        "One of the most common usage scenarios is to calculate the density of various",
621
        "groups across a lipid bilayer, typically with the z axis being the normal",
622
        "direction. For short simulations, small systems, and fixed box sizes this",
623
        "will work fine, but for the more general case lipid bilayers can be complicated.",
624
        "The first problem that while both proteins and lipids have low volume",
625
        "compressibility, lipids have quite high area compressiblity. This means the",
626
        "shape of the box (thickness and area/lipid) will fluctuate substantially even",
627
        "for a fully relaxed system. Since GROMACS places the box between the origin",
628
        "and positive coordinates, this in turn means that a bilayer centered in the",
629
        "box will move a bit up/down due to these fluctuations, and smear out your",
630
        "profile. The easiest way to fix this (if you want pressure coupling) is",
631
        "to use the [TT]-center[tt] option that calculates the density profile with",
632
        "respect to the center of the box. Note that you can still center on the",
633
        "bilayer part even if you have a complex non-symmetric system with a bilayer",
634
        "and, say, membrane proteins - then our output will simply have more values",
635
        "on one side of the (center) origin reference.[PAR]",
636

    
637
        "Even the centered calculation will lead to some smearing out the output",
638
        "profiles, as lipids themselves are compressed and expanded. In most cases",
639
        "you probably want this (since it corresponds to macroscopic experiments),",
640
        "but if you want to look at molecular details you can use the [TT]-relative[tt]",
641
        "option to attempt to remove even more of the effects of volume fluctuations.[PAR]",
642

    
643
        "Finally, large bilayers that are not subject to a surface tension will exhibit",
644
        "undulatory fluctuations, where there are 'waves' forming in the system.",
645
        "This is a fundamental property of the biological system, and if you are",
646
        "comparing against experiments you likely want to include the undulation",
647
        "smearing effect.",
648
        "",
649
    };
650

    
651
    output_env_t       oenv;
652
    static const char *dens_opt[] =
653
    { NULL, "mass", "number", "charge", "electron", NULL };
654
    static int         axis        = 2;  /* normal to memb. default z  */
655
    static const char *axtitle     = "Z";
656
    static int         nslices     = 50; /* nr of slices defined       */
657
    static int         ngrps       = 1;  /* nr. of groups              */
658
    static gmx_bool    bSymmetrize = FALSE;
659
    static gmx_bool    bCenter     = FALSE;
660
    static gmx_bool    bRelative   = FALSE;
661

    
662
    t_pargs            pa[]        = {
663
        { "-d", FALSE, etSTR, {&axtitle},
664
          "Take the normal on the membrane in direction X, Y or Z." },
665
        { "-sl",  FALSE, etINT, {&nslices},
666
          "Divide the box in this number of slices." },
667
        { "-dens",    FALSE, etENUM, {dens_opt},
668
          "Density"},
669
        { "-ng",       FALSE, etINT, {&ngrps},
670
          "Number of groups of which to compute densities." },
671
        { "-center",   FALSE, etBOOL, {&bCenter},
672
          "Perform the binning relative to the center of the (changing) box. Useful for bilayers." },
673
        { "-symm",     FALSE, etBOOL, {&bSymmetrize},
674
          "Symmetrize the density along the axis, with respect to the center. Useful for bilayers." },
675
        { "-relative", FALSE, etBOOL, {&bRelative},
676
          "Use relative coordinates for changing boxes and scale output by average dimensions." }
677
    };
678

    
679
    const char        *bugs[] = {
680
        "When calculating electron densities, atomnames are used instead of types. This is bad.",
681
    };
682

    
683
    double           **density;        /* density per slice          */
684
    real               slWidth;        /* width of one slice         */
685
    char              *grpname_center; /* centering group name     */
686
    char             **grpname;        /* groupnames                 */
687
    int                nr_electrons;   /* nr. electrons              */
688
    int                ncenter;        /* size of centering group    */
689
    int               *ngx;            /* sizes of groups            */
690
    t_electron        *el_tab;         /* tabel with nr. of electrons*/
691
    t_topology        *top;            /* topology               */
692
    int                ePBC;
693
    atom_id           *index_center;   /* index for centering group  */
694
    atom_id          **index;          /* indices for all groups     */
695
    int                i;
696

    
697
    t_filenm           fnm[] = { /* files for g_density       */
698
        { efTRX, "-f", NULL,  ffREAD },
699
        { efNDX, NULL, NULL,  ffOPTRD },
700
        { efTPR, NULL, NULL,  ffREAD },
701
        { efDAT, "-ei", "electrons", ffOPTRD }, /* file with nr. of electrons */
702
        { efXVG, "-o", "density", ffWRITE },
703
    };
704

    
705
#define NFILE asize(fnm)
706

    
707
    if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
708
                           NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs,
709
                           &oenv))
710
    {
711
        return 0;
712
    }
713

    
714
    if (bSymmetrize && !bCenter)
715
    {
716
        fprintf(stderr, "Can not symmetrize without centering. Turning on -center\n");
717
        bCenter = TRUE;
718
    }
719
    /* Calculate axis */
720
    axis = toupper(axtitle[0]) - 'X';
721

    
722
    top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC); /* read topology file */
723

    
724
    snew(grpname, ngrps);
725
    snew(index, ngrps);
726
    snew(ngx, ngrps);
727

    
728
    if (bCenter)
729
    {
730
        fprintf(stderr,
731
                "\nNote: that the center of mass is calculated inside the box without applying\n"
732
                "any special periodicity. If necessary, it is your responsibility to first use\n"
733
                "trjconv to make sure atoms in this group are placed in the right periodicity.\n\n"
734
                "Select the group to center density profiles around:\n");
735
        get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &ncenter,
736
                  &index_center, &grpname_center);
737
    }
738
    else
739
    {
740
        ncenter = 0;
741
    }
742

    
743
    fprintf(stderr, "\nSelect %d group%s to calculate density for:\n", ngrps, (ngrps > 1) ? "s" : "");
744
    get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, ngx, index, grpname);
745

    
746
    if (dens_opt[0][0] == 'e')
747
    {
748
        nr_electrons =  get_electrons(&el_tab, ftp2fn(efDAT, NFILE, fnm));
749
        fprintf(stderr, "Read %d atomtypes from datafile\n", nr_electrons);
750

    
751
        calc_electron_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density,
752
                              &nslices, top, ePBC, axis, ngrps, &slWidth, el_tab,
753
                              nr_electrons, bCenter, index_center, ncenter,
754
                              bRelative, oenv);
755
    }
756
    else
757
    {
758
        calc_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density, &nslices, top,
759
                     ePBC, axis, ngrps, &slWidth, bCenter, index_center, ncenter,
760
                     bRelative, oenv, dens_opt);
761
    }
762

    
763
    plot_density(density, opt2fn("-o", NFILE, fnm),
764
                 nslices, ngrps, grpname, slWidth, dens_opt,
765
                 bCenter, bRelative, bSymmetrize, oenv);
766

    
767
    do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");  /* view xvgr file */
768
    return 0;
769
}