Giter Club home page Giter Club logo

alchemical-analysis's Introduction

Alchemical Analysis: An open tool implementing some recommended practices for analyzing alchemical free energy calculations

Use Alchemlyb instead

We are in the process of migrating all functionality from here to instead use alchemlyb, which focuses on being a general Python library for performing analysis of alchemical calculations rather than a stand-alone command-line tool for analyzing these calculations. We recommend you move to using, and contribute to the development of, alchemlyb instead. It already has the following advantages, among others:

  • Python 3 compatible
  • Works as a library rather than a stand-alone tool, allowing easier connection with trajectory analysis
  • Has automated testing/continuous integration testing
  • Easily extensible

However, some of the plotting functionality and analysis features available in alchemical analysis are not yet available in alchemlyb. If there is functionality which is particularly important to you which is missing there, please raise an issue on the alchemlyb issue tracker.

In the meantime, Alchemical Analysis is still available but deprecated and we discourage its use. Please use alchemlyb instead. We may occasionally perform minor/critical maintenance here, but intend to cease maintaining this package over time.

About Alchemical Analysis

Analyze alchemical free energy calculations conducted in GROMACS, AMBER or SIRE using recommended best practices from Klimovich et al., JCAMD 29:397-411 (2015).

This tool handles analysis via a slate of free energy methods, including BAR, MBAR, TI, and the Zwanzig relationship (exponential averaging) among others, and provides a good deal of analysis of computed free energies and convergence in order to help you assess the quality of your results.

If you have problems with this tool, please use the github issue tracker.

Citation DOI for Citing Alchemical Analysis

Alchemical Analysis is research software. If you make use of it in work which you publish, please cite it. The BibTex reference is

@article{Klimovich:2015er,
author = {Klimovich, Pavel V and Shirts, Michael R and Mobley, David L},
title = {Guidelines for the analysis of free energy calculations},
journal = {J Comput Aided Mol Des},
year = {2015},
volume = {29},
number = {5},
pages = {397--411},
doi = {10.1007/s10822-015-9840-9}
}

Prerequisites

Alchemical Analysis requires the pymbar module (we recommend installation via conda, but other options are also available). Version 3.0.0.dev0 or above is required. The latest version of pymbar can be found here.

Python 2 is required.

Installation

git clone https://github.com/MobleyLab/alchemical-analysis.git
cd alchemical-analysis
sudo python setup.py install

Usage

Script: alchemical_analysis

This implements recommended practices for analyzing alchemical free energy calculations, as described in Klimovich et al., JCAMD 29:397-411 (2015). This was motivated in part by earlier work illustrating how to apply MBAR to alchemical free energy calculations (and a comparison with other methods) in Paliwal and Shirts, J. Chem. Theory Comp, v. 7, 4115-4134 (2011).

See description in samples/gromacs/README.md, samples/sire/README.md, and samples/amber/README.md.

Help for alchemical_analysis.py (obtained with alchemical_analysis -h) is:

  -h, --help            show this help message and exit
  -a SOFTWARE, --software=SOFTWARE
                        Package's name the data files come from: Gromacs,
                        Sire, Desmond, or AMBER. Default: Gromacs.
  -b FRACTION, --fraction=FRACTION
                        The fraction of the energy file will be used to
                        calculate the statistics. Default: 1.0
  -c, --cfm             The Curve-Fitting-Method-based consistency inspector.
                        Default: False.
  -d DATAFILE_DIRECTORY, --dir=DATAFILE_DIRECTORY
                        Directory in which data files are stored. Default:
                        Current directory.
  -e, --backward        Extract the energy data from the backward direction.
                        Default: False
  -f BFORWREV, --forwrev=BFORWREV
                        Plot the free energy change as a function of time in
                        both directions, with the specified number of points
                        in the time plot. The number of time points (an
                        integer) must be provided. Default: 0
  -g, --breakdown       Plot the free energy differences evaluated for each
                        pair of adjacent states for all methods, including the
                        dH/dlambda curve for TI. Default: False.
  -i UNCORR_THRESHOLD, --threshold=UNCORR_THRESHOLD
                        Proceed with correlated samples if the number of
                        uncorrelated samples is found to be less than this
                        number. If 0 is given, the time series analysis will
                        not be performed at all. Default: 50.
  -j RESULTFILENAME, --resultfilename=RESULTFILENAME
                        custom defined result filename prefix. Default:
                        results
  -k BSKIPLAMBDAINDEX, --koff=BSKIPLAMBDAINDEX
                        Give a string of lambda indices separated by '-' and
                        they will be removed from the analysis. (Another
                        approach is to have only the files of interest present
                        in the directory). Default: None.
  -m METHODS, --methods=METHODS
                        A list of the methods to esitimate the free energy
                        with. Default: [TI, TI-CUBIC, DEXP, IEXP, BAR, MBAR].
                        To add/remove methods to the above list provide a
                        string formed of the method strings preceded with +/-.
                        For example, '-ti_cubic+gdel' will turn methods into
                        [TI, DEXP, IEXP, BAR, MBAR, GDEL]. 'ti_cubic+gdel', on
                        the other hand, will call [TI-CUBIC, GDEL]. 'all'
                        calls the full list of supported methods [TI, TI-
                        CUBIC, DEXP, IEXP, GINS, GDEL, BAR, UBAR, RBAR, MBAR].
  -n UNCORR, --uncorr=UNCORR
                        The observable to be used for the autocorrelation
                        analysis; either 'dhdl_all' (obtained as a sum over
                        all energy components) or 'dhdl' (obtained as a sum
                        over those energy components that are changing;
                        default) or 'dE'. In the latter case the energy
                        differences dE_{i,i+1} (dE_{i,i-1} for the last
                        lambda) are used.
  -o OUTPUT_DIRECTORY, --out=OUTPUT_DIRECTORY
                        Directory in which the output files produced by this
                        script will be stored. Default: Same as
                        datafile_directory.
  -p PREFIX, --prefix=PREFIX
                        Prefix for datafile sets, i.e.'dhdl' (default).
  -q SUFFIX, --suffix=SUFFIX
                        Suffix for datafile sets, i.e. 'xvg' (default).
  -r DECIMAL, --decimal=DECIMAL
                        The number of decimal places the free energies are to
                        be reported with. No worries, this is for the text
                        output only; the full-precision data will be stored in
                        'results.pickle'. Default: 3.
  -s EQUILTIME, --skiptime=EQUILTIME
                        Discard data prior to this specified time as
                        'equilibration' data. Units picoseconds. Default: 0
                        ps.
  -t TEMPERATURE, --temperature=TEMPERATURE
                        Temperature in K. Default: 298 K.
  -u UNITS, --units=UNITS
                        Units to report energies: 'kJ', 'kcal', and 'kBT'.
                        Default: 'kJ'
  -v, --verbose         Verbose option. Default: False.
  -w, --overlap         Print out and plot the overlap matrix. Default: False.
  -x, --ignoreWL        Do not check whether the WL weights are equilibrated.
                        No log file needed as an accompanying input.
  -y RELATIVE_TOLERANCE, --tolerance=RELATIVE_TOLERANCE
                        Convergence criterion for the energy estimates with
                        BAR and MBAR. Default: 1e-10.
  -z INIT_WITH, --initialize=INIT_WITH
                        The initial MBAR free energy guess; either 'BAR' or
                        'zeros'. Default: 'BAR'.

License

MIT.

alchemical-analysis's People

Contributors

davidlmobley avatar fepanalysis avatar halx avatar harlor avatar kyleabeauchamp avatar mrshirts avatar shuail avatar wesbarnett avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

alchemical-analysis's Issues

Fix AMBER parser to use units consistent with the rest of the code

See #45 , @halx : The output units are inconsistent. For example, output shows:

The gradients from the correlated samples (kcal/mol):
0.00000 -15.817
0.10000 -24.593
0.20000 -32.826
0.30000 -44.704
0.40000 -56.149
0.50000 -61.926
0.60000 -65.991
0.70000 -77.821
0.80000 -103.103
0.90000 -134.217
1.00000 -125.538
TI = -67.201

but later

Estimating the free energy change with TI...

States TI (kJ/mol)

0 -- 1 -8.454 +- 0.333 
1 -- 2 -12.068 +- 0.506 
2 -- 3 -16.160 +- 0.597 
3 -- 4 -20.983 +- 0.569 
4 -- 5 -24.701 +- 0.445 
5 -- 6 -26.775 +- 0.432 
6 -- 7 -30.216 +- 0.564 
7 -- 8 -37.965 +- 0.551 
8 -- 9 -49.648 +- 0.448 
9 -- 10 -54.341 +- 0.337

TOTAL:     -281.312  +-  2.176 

Using two different sets of units is confusing and can easily cause the user to think the results are inconsistent.

Corrupt xvgs check slows down gromacs parser and makes unnecessary backups

The addition of checking for corrupt xvgs slows down the gromacs parser quite a bit which may or may not be related to the use of numpy.genfromtxt? It might be better/faster to iterate over the lines and check for lines with different numbers of elements.

In the gromacs parser, prior to the call removeCorruptLines(), I had it make a backup of the xvgs in a separate folder. This is not always totally necessary in the case where there are no corrupt lines in the xvgs. Instead, it should only make backups of offending xvgs.

multiple expanded ensemble files not working

One feature the code use to have was the ability to load together multiple expanded ensemble files. For example, there might be 18 states, but 4 separate expanded ensemble replicas. Previously (several generations ago) the code could load in all 4 of those files, and put all the corresponding samples together into a single u_kln matrix. However, after the edits, this can't be done.

One issue that's relatively easily corrected is for expanded ensemble reading the number of states from the file, rather than counting the number of files (which may not be the same). I've posted a pull_request with this change. However, it still is making some assumptions about how the data is loaded in to the dhdl array and u_kln matrix, and I don't know that I can updated those without messing up the data processing changes.

Any thoughts for how this functionality could be readded?

Error in Reading GROMACS dhdl.xvg with NumPy v1.10

Hello:

I have used alchemical_analysis.py in the past successfully to estimate MBAR free energy differences between states. Thank you for making this tool available for open use.

I recently updated NumPy to version 1.10. In trying to use alchemical_analysis.py on a GROMACS dhdl.xvg, I received the following error:

File "~/git/alchemical-analysis/alchemical_analysis/parser_gromacs.py", line 302, in readDataGromacs
extract_states = numpy.genfromtxt(f.filename, dtype=int, skiprows=f.skip_lines, skip_footer=1*bLenConsistency, usecols=1)
TypeError: genfromtxt() got an unexpected keyword argument 'skiprows'

Looking at the NumPy documentation at http://docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.genfromtxt.html,

numpy.genfromtxt

numpy.genfromtxt(...)
Load data from a text file, with missing values handled as specified.

Parameters:
skiprows: int, optional
skiprows was removed in numpy 1.10. Please use skip_header instead.
skip_header: int, optional
The number of lines to skip at the beginning of the file.

In parser_gromacs.py, the following change of line 302 fixes the present issue:
extract_states = numpy.genfromtxt(f.filename, dtype=int, skip_header=f.skip_lines, skip_footer=1*bLenConsistency, usecols=1)

This may not be all the changes due to NumPy v1.10, but it does in this case allow the code to progress farther. I run into another issue, which may not be related to the NumPy update.
After the change for a single .xvg, I incur the following error:

File "~/git/alchemical-analysis/alchemical_analysis/parser_gromacs.py", line 126, in slice_data
dhdlt[state, :, nsnapshots_l[state]:nsnapshots_r[state]] = data[read_dhdl_sta : read_dhdl_end, :]
ValueError: could not broadcast input array from shape (3,1808) into shape (3,10001)

I'm not entirely sure why this error occurs. I've only seen broadcasting errors like this in trying to run multiple .xvgs which are not each from a unique state (Related to now-closed issue multiple expanded ensemble files not working #36). There are 10001 data points in the xvg, and I presume the 3 is from the three dhdl potential types fep, coul, and vdw. 1808 appears to be how many samples are in state 0. Just listing the details here in case this error is related to the NumPy version.

Each of these errors (genfromtxt TypeError and dhdlt ValueError) were repeated after pulling the latest version of the master branch (immediately before posting this issue). The specific NumPy version is 1.10.4-py27_0, installed by conda.

I'm not sure if there is an intent to keep up with newer versions of NumPy, but I wanted to make a note of the error, in case it helps.

Thank you for your time.

Issue with Sorting input data by starting time (parser_amber.py)

Hello there,

I am getting this error when I try to use alchemical_analysis.py which is coming from amber MDOUT file parsing. I had it write out the offending variables ( found, uniq in parser_amber.py)

Command line was: ./alchemical-analysis/alchemical_analysis/alchemical_analysis.py -a AMBER -d data -p [01].*/ti00[2-9] -q out -o . -r 5 -u kcal -s 0 -g -w

Loading in data from data/0.4/ti003.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/0.4/ti002.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/0.7/ti003.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/0.7/ti002.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/0.6/ti003.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/0.6/ti002.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/0.0/ti003.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/0.0/ti002.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/0.9/ti003.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/0.9/ti002.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/0.2/ti003.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/0.2/ti002.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/0.1/ti003.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/0.1/ti002.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/0.5/ti003.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/0.5/ti002.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/0.3/ti003.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/0.3/ti002.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/0.8/ti003.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/0.8/ti002.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/1.0/ti003.out... 6000 data points, 6000 DV/DL averages
Loading in data from data/1.0/ti002.out... 6000 data points, 6000 DV/DL averages
Sorting input data by starting time
found: [200.0, 200.0] 2 uniq: set([200.0]) 1
ERROR: Same starting time occurs multiple times.

The header for parser_amber.py says

"Relevant namelist variables in MDIN: ntpr, ntave, ifmbar and bar*,
do NOT set t as the code depends on a correctly set begin time"

and it appears the issue stems from the begin time mentioned here. I am not exactly sure what this time it mentions here is.

Any thoughts appreciated.

Thanks!

Emmanuel

Inaccurate decorrelation times when lambda_coul and lambda_vdw change at different times.

We performed a calculation turning off an ion with separate lambda_coul and lambda_vdw, and where we first decoupled the coulomb and then the LJ. When alchemical-analysis calculated uncorrelated data points, some of the lambda_coul = 0 / lambda_LJ changing states appeared to have very long correlation times. However, it was computing the correlation time using the sum of dH/dlambda_coul and dH/dlambda_vdw at all lambdas.

Because the LJ soft core was coming off, the lambda_coul spiked quite a bit at an infrequent interval, leading to long correlation times. However, this is spurious, since the dH/dlambda_coul is not used in the calculation of the free energy change.

Upon further inspection, the fact that the sum of the two dH/dlambdas is used affects the lambda_coul = 0 / lambda_LJ = 0 end as well since the decorrelation time is dominated by very fast fluctuations in dH/dlambda_coul, whereas dH/dlambda_LJ has slower fluctuations (which are being washed out).

One possible solution is to only sum together the decorrelation times the lambdas that are changing, since delta G = delta lambda \dot gradient H. That seems a better solution (it matches the observable, where ), but I wanted to verify that other developers agree before trying to make the change. Any other ideas?

General suggestions for code cleanup/improvements

There is currently a lot of cruft growing in the code which is basically down to design choices.

Here some suggestions:
The parsers should be designed in a "polymorphic" fashion such that the main code does not need to know anything about the individual parser codes. This implies that a clear interface and a specification must be designed to define what data is passed in and what is passed out and in what form. It is the responsibility of the parser code to signal back what it can and cannot do. E.g., as it stands now, dhdlt should be None if MD code can't do TI or u_klt should be None if it *BAR analysis can't be done. Energies should be converted to a "canonical" form or it should be signalled back to the main code what the units are. This will relieve the main code maintainer to do code-specific clean-up operations.

Optparse should be replaced with argparse.

The globals need to go as this is really bad programming style. The problem with this is that all these variables are accessible throughout all of the module which mean that they can be modified anywhere. This can create hard to track down bugs (all you need is a second person messing around there...).

order of methods

I have just noticed that the order of methods is
..., 'DEXP', 'IEXP', 'GINS', 'GDEL', ...

I guess the D stands for DEL and I for INS. So the order is exactly the opposite between the EXP formula and the Gaussian approximation. Is this on purpose?

I also wonder what INSert and DELete stand for? Does this terminology follow from absolute transformations? In relative ones you can have both.

-a amber -m all throws error

Hi David,

I tried alchemical_analysis.py with the amber test files provided in samples/amber/data

with
python2.7 ../../../alchemical_analysis.py -a amber -p "ti*/ti" -q out -o ~/test

the script runs through and produces the expected output for TI and TI-CUBIC.

if I add the option -m all to the above call, I get the following traceback:

Traceback (most recent call last):
File "../../../alchemical_analysis/alchemical_analysis.py", line 1184, in main()
File "../../../alchemical_analysis/alchemical_analysis.py", line 1160, in main
Deltaf_ij, dDeltaf_ij = estimatewithMBAR(u_kln, N_k, P.relative_tolerance, regular_estimate=True)
NameError: global name 'u_kln' is not defined

Is this a bug or did I do something wrong? Any help appreciated.

Thanks,
Th.

index out of bounds error with some values of bForwrev

I find that when using certain values with the flag "-f", I get an index out of bounds error because the array ts in the function dF_t() is one element longer than expected. I found this issue with -f values of 6 or 28 so far.

I believe it has to do with the fact the array tf is defined with numpy.arange:
tf = numpy.arange(0,1+increment,increment)(max(nsnapshots)-1)+1
In some cases (like for 6 or 28) it seems that the 'stop' (1+increment) element is included in the array, later on resulting in a ts array that is one element longer than expected, and with ts[-1] having a time over the end time of the simulation. I imagine it could be a problem related to the representation of floating numbers, and possibly replacing arange with linspace might be enough to solve the issue:
tf = numpy.linspace(0.0, 1.0, n_tf)
(max(nsnapshots)-1)+1

error with -c -m ti+ti-cubic

The traceback below suggest to me that the code doesn't properly guard against this combination of command line switches. Obviously it doesn't make sense to do the CFM analysis when there is not BAR type data.

Traceback (most recent call last):
File "/home/hhl/projects/alchemical-analysis/alchemical_analysis/alchemical_analysis.py", line 1224, in
main()
File "/home/hhl/projects/alchemical-analysis/alchemical_analysis/alchemical_analysis.py", line 1219, in main
plotCFM(u_kln, N_k, 50)
File "/home/hhl/projects/alchemical-analysis/alchemical_analysis/alchemical_analysis.py", line 1111, in plotCFM
plotdg_vs_dU(yy, df_allk, ddf_allk)
File "/home/hhl/projects/alchemical-analysis/alchemical_analysis/alchemical_analysis.py", line 1046, in plotdg_vs_dU
ax.fill_between(xx_i, df_allk[i]['BAR'] - ddf_allk[i]['BAR'], df_allk[i]['BAR'] + ddf_allk[i]['BAR'], color='#D2B9D3', zorder=-1)
KeyError: 'BAR'

index 40 is out of bounds for axis 0 with size 40

I am calculating binding free energy of Protein ligand complex, 41 lambdas. i got this error.
Traceback (most recent call last):
File "/share/anaconda2/bin/alchemical_analysis.py", line 4, in
import('pkg_resources').run_script('alchemical-analysis==1.0.2.dev0', 'alchemical_analysis.py')
File "/share/anaconda2/lib/python2.7/site-packages/setuptools-20.3-py2.7.egg/pkg_resources/init.py", line 726, in run_script
File "/share/anaconda2/lib/python2.7/site-packages/setuptools-20.3-py2.7.egg/pkg_resources/init.py", line 1484, in run_script
File "/share/anaconda2/lib/python2.7/site-packages/alchemical_analysis-1.0.2.dev0-py2.7.egg/EGG-INFO/scripts/alchemical_analysis.py", line 1226, in
main()
File "/share/anaconda2/lib/python2.7/site-packages/alchemical_analysis-1.0.2.dev0-py2.7.egg/EGG-INFO/scripts/alchemical_analysis.py", line 1212, in main
totalEnergies()
File "/share/anaconda2/lib/python2.7/site-packages/alchemical_analysis-1.0.2.dev0-py2.7.egg/EGG-INFO/scripts/alchemical_analysis.py", line 583, in totalEnergies
dF[name] += df_allk[k][name]
IndexError: index 40 is out of bounds for axis 0 with size 40

Consider implementing extrapolation scheme for MBAR to handle AMBER's lack of endpoint data

I wonder if it would make sense to introduce an extrapolation scheme for "incomplete" MBAR data data. The principal problem for me is sander which cannot sample the end-points. I also heard from some other person that he doesn't sample the end-points and relies on extrapolation apparently. Is there any good reason why that should not be done? I certainly understand that this shouldn't be used to bridge larger gaps.

gromacs parser

i am a beginner, can you guide me about this error, how should i fix it.

Traceback (most recent call last):
File "alchemical_analysis.py", line 1076, in
import parsers.parser_gromacs
ImportError: No module named parsers.parser_gromacs

Generalize handling of other (non-charging/vdW) free energy components, such as for restraints, so that they don't get grouped with vdW

A possible small enhancement: since for absolute calculations there is a also a restraints component other than coulomb and vdw, it could be good to have the free energy of this term shown as well in the summary table of the results. At the moment I think the restraints free energy ends up being included in the vdw component, as suggested by the remark.

It might be possible to do this just by modifying slightly the totalEnergies() function? I attach a snippet of quickly modified code which seemed to work for me; it could have some issues though, so it is really just an idea.
showrestraints.txt

Sample SIRE data doesn't process properly

I ran the latest master alchemical_analysis on the sire data under samples as per the relevant README, and ran into the below error. @halx , can you have a look?

Loading in data from data/actual_grad_0000.dat  (state 0) ...
Traceback (most recent call last):
  File "../../alchemical_analysis/alchemical_analysis.py", line 1224, in <module>
    main()
  File "../../alchemical_analysis/alchemical_analysis.py", line 1169, in main
    nsnapshots, lv, dhdlt, u_klt = parser_sire.readDataSire(P)
  File "/Users/dmobley/GoogleDrive/research/github/alchemical-analysis/alchemical_analysis/parser_sire.py", line 90, in readDataSire
    f.loadtxtSire(nf)
  File "/Users/dmobley/GoogleDrive/research/github/alchemical-analysis/alchemical_analysis/parser_sire.py", line 63, in loadtxtSire
    dhdlt[state, :, :nsnapshots[state]] = numpy.genfromtxt(self.filename, dtype=float, skiprows=self.skip_lines, usecols=1)
TypeError: genfromtxt() got an unexpected keyword argument 'skiprows'

Actually, a quick Google suggests that this is a problem with an updated numpy. I may be able to quickly fix. Specifically:

skiprows : int, optional
skiprows was removed in numpy 1.10. Please use skip_header instead.

Lambda index out of bounds when accumulating free energies

Hi,

I'm getting the following error when I use alchemical-analysis for the second part of my run (bonded, LJ and masses).

File "/Users/dvr/Work/alchemical-analysis-master/alchemical_analysis/alchemical_analysis.py", line 1226, in
main()
File "/Users/dvr/Work/alchemical-analysis-master/alchemical_analysis/alchemical_analysis.py", line 1212, in main
totalEnergies()
File "/Users/dvr/Work/alchemical-analysis-master/alchemical_analysis/alchemical_analysis.py", line 583, in totalEnergies
dF[name] += df_allk[k][name]
IndexError: index out of bounds

Everything works fine for the first and third part of the transformation with the same alchemical-analysis version and settings. Detailed output and input files can be found here:
https://www.dropbox.com/sh/vmjh7rmswbg0bd8/AADUpt8622PBUnxill6eOoIma?dl=0

Kind regards

Dries

Make an official (numbered) release version of the code published with our paper, then begin following version numbering conventions

See http://semver.org/ for info on version numbering conventions.

We should make an official (numbered) release version corresponding to the code published with our paper, and clearly indicate in the version description what this is. We then ought to begin tracking "official" versions so these can be specified by people writing papers using this code. Issues can then be tied to specific versions for their resolution.

MBAR method fails

Hi there,

I was trying to calculate the binding affinity by using the MBAR method but unfortunately it fails with the error:

Traceback (most recent call last):
File "/home/gac/local/Common/bin/alchemical_analysis.py", line 1104, in
Deltaf_ij, dDeltaf_ij = estimatewithMBAR(u_kln, N_k, P.relative_tolerance, regular_estimate=True)
File "/home/gac/local/Common/bin/alchemical_analysis.py", line 240, in estimatewithMBAR
(Deltaf_ij, dDeltaf_ij, theta_ij) = MBAR.getFreeEnergyDifferences(uncertainty_method='svd-ew')
ValueError: need more than 2 values to unpack

I don't have any problem by selecting TI or BAR.

The link related to the directory that I would like to process is:

https://drive.google.com/open?id=0B4Hyb-RxLkeGTEl4QVltT3NhYkk

Cheers,

Gaetano

Coulomb Free Energy

The Coulomb part of the free energy calculation is being reported as zero for hydration and solvation free energy calculations when we wouldn't expect it to be. This doesn't seem to be happening with the sample files that are provided. I can provide a sample set of .xvg files if needed.

Overlap matrix calculation doesn't work with current pymbar

Output of the overlap matrix seems to fail with current (2.1 and later - i.e. May 2014) versions of pymbar. @FEPanalysis , can you take a look? Apparently the default output of MBAR.computeOverlap() is now a scalar. Superficially it looks like switching the output to 'matrix' ought to do the trick, i.e. O = MBAR.computeOverlap(output='matrix')[2]', though this still seems to fail because of something deeper. I'm guessing the doc string of computeOverlap will tell you what needs to be changed, but I don't apparently have the version of pymbar you were using for this so it's unclear to me how the returns have changed.

NameError: global name 'parametersB' is not defined

I was using a different format of topologies files and I got this error saying the variable parametersB is not defined previously. The issue comes from Part IV, when it is processing the dihedrals and the top file format doesn't match with what the code is expecting, going to an exception in line 380. What I have noticed is acpype and parmed return different dihedrals functions types (9 and 1), and amb2gmx.pl returns the function type 3 (the one that works). I can provide samples in case it is necessary.

Traceback (most recent call last):
File "../../local/alchemical-setup-master/alchemical-setup.py", line 549, in
buildBonded() # Part IV.
File "../../local/alchemical-setup-master/alchemical-setup.py", line 429, in buildBonded
postprocessDihedrals()
File "../../local/alchemical-setup-master/alchemical-setup.py", line 385, in postprocessDihedrals
if l['parameters'][-1] != parametersB[-1]:
NameError: global name 'parametersB' is not defined

Provide a utility tool to attempt cleanup of bad GROMACS data

On rare occasions, GROMACS writes get interrupted by networking issues resulting on corrupted data on one line of an otherwise complete xvg file used by alchemical analysis. Currently this results in an exception (seen below). Reworking the code to deal with this would involve major efficiency hits as far as I can tell. Probably a better solution, at least for now, is to provide a utility tool which can be run to clean up the .xvg files in question in problematic cases (basically by deleting corrupted lines - which typically seem to be just a single line).

Here's an example set of data:
https://drive.google.com/file/d/0B5sFkoaaNUxtRWtYWTY3RnBUMnM/view?usp=sharing

Here's the exception:
Loading in data from ./prod.1.xvg (state 1) ...
Traceback (most recent call last):
File "/home/pklimovi/analysis/alchemical-free-energy/alchemical_analysis.py", line 1077, in
nsnapshots, lv, dhdlt, u_klt = parsers.parser_gromacs.readDataGromacs(P)
File "/data11/home/pklimovi/analysis/alchemical-free-energy/parsers/parser_gromacs.py", line 280, in readDataGromacs
f.iter_loadtxt(nf)
File "/data11/home/pklimovi/analysis/alchemical-free-energy/parsers/parser_gromacs.py", line 120, in iter_loadtxt
data = numpy.fromiter(iter_func(), dtype=float)
ValueError: invalid literal for float(): 0.8063

If I search for that string in the file it looks like there's a bad (corrupted) line in the file on line 12520.

Properly handle case of all zeros in all input .xvg files

Currently in this case the script fails with an exception due to a problem calculating the statistical inefficiency:

Traceback (most recent call last):
File "/home/pklimovi/analysis/alchemical-free-energy/alchemical_analysis.py", line 1099, in
dhdl, N_k, u_kln = uncorrelate(sta=numpy.zeros(K, int), fin=nsnapshots, do_dhdl=True)
File "/home/pklimovi/analysis/alchemical-free-energy/alchemical_analysis.py", line 160, in uncorrelate
g[k] = pymbar.timeseries.statisticalInefficiency(dhdl_sum)
File "/home/pklimovi/analysis/pymbar-2.1.0-beta/pymbar/timeseries.py", line 178, in statisticalInefficiency
raise ParameterError('Sample covariance sigma_AB^2 = 0 -- cannot compute statistical inefficiency')
pymbar.utils.ParameterError: Sample covariance sigma_AB^2 = 0 -- cannot compute statistical inefficiency

See for example /work/cluster/vjaber/normal_hydrations/CRASHED_TRAJ/1234_tetrachlorobenzene_to_4_chlorophenol_vac

We ought instead to catch this case earlier (i.e. when the data is loaded from the .xvg files) and print out a note along with a sensible results.txt file. (This is important in, for example, @Limn1 's calculations since some calculations which yield zeros were accidentally run.)

Problem with the Curve-Fitting-Method-based consistency inspector

As suggested I move my post from Pymbar to here.

I am running pymbar 3.0 with the alchemical_analysis without problems except for the " -c " flag.

On the provided sample input data it is running fine but if I want to use it on my dhdl files I get the following mistake:

Plotting the CFM figure...
/opt/soft/anaconda-2.2.0/lib/python2.7/site-packages/matplotlib/axes/_base.py:2562: UserWarning: Attempting to set identical left==right results
in singular transformations; automatically expanding.
left=17.7338848684, right=17.7338848684
'left=%s, right=%s') % (left, right))
/opt/soft/anaconda-2.2.0/lib/python2.7/site-packages/matplotlib/axes/_base.py:2562: UserWarning: Attempting to set identical left==right results
in singular transformations; automatically expanding.
left=17.7624749232, right=17.7624749232
'left=%s, right=%s') % (left, right))

In case this is of interest: The dhdl-files were generated with Gromacs 5.0.4

My files are at

https://drive.google.com/folderview?id=0B2M9aqeJrxnYYTlsdU9jMnBQOGs&usp=sharing

Best regards,
Julian

cannot use MBAR

Dear Prof Mobley
"Number of states is NOT the same for all simulations; I'm assuming that we only evaluate
nearest neighbor states, and so cannot use MBAR, removing the method."

does it mean that i did mistake in numbr\ering states in mdp files? all the mdp files are same except the init_lambda_state values.

Change README.md to properly point

Currently the README.md refers to:

gromacs/README.md, sire/README.md, and amber/README.md

But these folders don't exist. I think it should be

samples/gromacs/README.md

etc.

small ligand solvation energy calculation:error:cannot compute statistical inefficiency

i have done ligand simulation in water and now wantto calculate solvation free energy , i used 31 lambdas states, i got this error by using alchemical_analysis.py.

Traceback (most recent call last):
File "/share/anaconda2/bin/alchemical_analysis.py", line 4, in
import('pkg_resources').run_script('alchemical-analysis==1.0.2.dev0', 'alchemical_analysis.py')
File "/share/anaconda2/lib/python2.7/site-packages/setuptools-19.6.2-py2.7.egg/pkg_resources/init.py", line 724, in run_script
File "/share/anaconda2/lib/python2.7/site-packages/setuptools-19.6.2-py2.7.egg/pkg_resources/init.py", line 1649, in run_script
File "/share/anaconda2/lib/python2.7/site-packages/alchemical_analysis-1.0.2.dev0-py2.7.egg/EGG-INFO/scripts/alchemical_analysis.py", line 1226, in
main()
File "/share/anaconda2/lib/python2.7/site-packages/alchemical_analysis-1.0.2.dev0-py2.7.egg/EGG-INFO/scripts/alchemical_analysis.py", line 1197, in main
dhdl, N_k, u_kln = uncorrelate(sta=numpy.zeros(K, int), fin=nsnapshots, do_dhdl=True)
File "/share/anaconda2/lib/python2.7/site-packages/alchemical_analysis-1.0.2.dev0-py2.7.egg/EGG-INFO/scripts/alchemical_analysis.py", line 181, in uncorrelate
g[k] = pymbar.timeseries.statisticalInefficiency(dhdl_sum)
File "build/bdist.linux-x86_64/egg/pymbar/timeseries.py", line 178, in statisticalInefficiency
pymbar.utils.ParameterError: Sample covariance sigma_AB^2 = 0 -- cannot compute statistical inefficiency

cannot use MBAR

Dear Prof Mobley
"Number of states is NOT the same for all simulations; I'm assuming that we only evaluate
nearest neighbor states, and so cannot use MBAR, removing the method."

does it mean that i did mistake in numbr\ering states in mdp files? all the mdp files are same except the init_lambda_state values.

Pymbar can't estimate bound on free energy

I have performed a bunch of thermodynamic-integration simulations for liquid water in Gromacs, and am trying to analyze the .xvg files. I'm using Windows 10 and installed Anaconda2 (64 bit) with the gui installer and then pymbar through Anaconda2 before downloading alchemical-analysis-master.

In a Windows command prompt, I entered:

python alchemical_analysis.py -d C:<xvg_directory>\ -t 200 -p 1 -v > results

This is the result:
Traceback (most recent call last):
File "alchemical_analysis.py", line 1224, in
main()
File "alchemical_analysis.py", line 1166, in main
nsnapshots, lv, dhdlt, u_klt = parser_gromacs.readDataGromacs(P)
File "C:<...>\alchemical-analysis-master\alchemical_analysis\parser_gromacs.py", line 187, in readDataGromacs
fs = sorted(fs, key=F.sortedHelper)
File "C:<...>\alchemical-analysis-master\alchemical_analysis\parser_gromacs.py", line 48, in sortedHelper
self.state = l[0] = int(l[0]) # Will be of use for selective MBAR analysis.
IndexError: list index out of range

where <...> is just the directory to which I downloaded alchemical-analysis-master. I don't understand the error and cannot find help for it online. In the directory where I have the .xvg files, all the files are simply named integers: 0.xvg, 1.xvg, 2.xvg, .... 20.xvg, and there are no other files in that folder. Am I calling the script correctly? Could it be the operating system (Windows vs. Linux)? I might be able to get on a Linux machine if it would help.

Long-term: Implement automated testing of alchemical-analysis

We should automatically check to ensure that it works as expected with new pymbar updates.

Use cases to check:

  • Expanded ensemble
  • Independent lambda values with and without charging/vdW transformations
    • with a single lambda vector
    • with multiple lambda vectors
  • Lambda exchange
  • With and without calc-lambda-neighbors

analysis TI with alchemical_analysis.py errors

Dear Prof. Mobley,

I run a TI analysis with AMBER and I installed and I am trying to use your alchemical_analysis.py script for the example that comes with the with it, from the sampes/amber directory as it is. The scripts runs, but I get an error I don’t understand, probably my python ignorance, here it is the ouptut, it seems it starts to run, but it gives a mbar_lambdas error.

/Users/admin/Data/TOOLS/alchemical-analysis-master/samples/amber
03:24 PM Sun Dec 13$ ../../alchemical_analysis/alchemical_analysis.py -a AMBER -d data/* -m TI -p ti00* -q out

Started on Sun Dec 13 15:25:08 2015

Command line was: ../../alchemical_analysis/alchemical_analysis.py -a AMBER -d data/0.0 data/0.1 data/0.2 data/0.3 data/0.4 data/0.5 data/0.6 data/0.7 data/0.8 data/0.9 data/1.0 -m TI -p ti00* -q out

Loading in data from data/0.0/ti002.out... 50 data points, 10 DV/DL averages
Loading in data from data/0.0/ti003.out... 50 data points, 10 DV/DL averages

Sorting input data by starting time and lambda
Traceback (most recent call last):
File "../../alchemical_analysis/alchemical_analysis.py", line 1224, in
main()
File "../../alchemical_analysis/alchemical_analysis.py", line 1172, in main
nsnapshots, lv, dhdlt, u_klt = parser_amber.readDataAmber(P)
File "/Users/admin/Data/TOOLS/alchemical-analysis-master/alchemical_analysis/parser_amber.py", line 582, in readDataAmber
if len(dvdl_all) != len(mbar_lambdas):
UnboundLocalError: local variable 'mbar_lambdas' referenced before assignment

If I use your script run.sh I get a different error:

./run.sh
Traceback (most recent call last):
File "../../alchemical_analysis//alchemical_analysis.py", line 1224, in
main()
File "../../alchemical_analysis//alchemical_analysis.py", line 1200, in main
Deltaf_ij, dDeltaf_ij = estimatewithMBAR(u_kln, N_k, P.relative_tolerance, regular_estimate=True)
File "../../alchemical_analysis//alchemical_analysis.py", line 290, in estimatewithMBAR
O = MBAR.computeOverlap()[2]
IndexError: invalid index to scalar variable.
/Users/admin/Data/TOOLS/alchemical-analysis-master/samples/amber
03:41 PM Sun Dec 13$ more run.sh
../../alchemical_analysis//alchemical_analysis.py
-a AMBER
-d data
-p '[01].*/ti00[2-9]'
-q out
-o .
-r 5
-u kcal
-s 0
-g
-w \

run.log
/Users/admin/Data/TOOLS/alchemical-analysis-master/samples/amber

Could you please help me understand what I am doing wrong?

Thanks a lot in advance,

Fabian

Dr. Fabian Glaser
Head of the Structural Bioinformatics section

Bioinformatics Knowledge Unit - BKU
The Lorry I. Lokey Interdisciplinary Center for Life Sciences and Engineering
Technion - Israel Institute of Technology, Haifa 32000, ISRAEL

fglaser at technion dot ac dot il
Tel: +972 4 8293701
http://bku.technion.ac.il

Problem with pymbar API; issues with cubic spline interpolation in TI

Hi,

I'm experiencing the following error when I try to calculate the overlap matrix:

File "/Users/dvr/Documents/Cheminformatics_tools/alchemical-analysis-master/alchemical_analysis/alchemical_analysis.py", line 1226, in
main()
File "/Users/dvr/Documents/Cheminformatics_tools/alchemical-analysis-master/alchemical_analysis/alchemical_analysis.py", line 1202, in main
Deltaf_ij, dDeltaf_ij = estimatewithMBAR(u_kln, N_k, P.relative_tolerance, regular_estimate=True)
File "/Users/dvr/Documents/Cheminformatics_tools/alchemical-analysis-master/alchemical_analysis/alchemical_analysis.py", line 292, in estimatewithMBAR
O = MBAR.computeOverlap()[2]
IndexError: invalid index to scalar variable.

My files can be found at https://www.dropbox.com/sh/sobp5423rm0ubog/AAAmrsxBLSXQVoein7ue9lZTa?dl=0. The command I'm using is: python ~/Documents/Cheminformatics_tools/alchemical-analysis-master/alchemical_analysis/alchemical_analysis.py -q xvg -p md -u kJ -c -t 300 -w.

This is only a prototype run to check if my lambda spacing is sufficient, hence the short simulations. The cfm.pdf looks weird as well, though I suppose that could be related to the length of the simulations. This is the first step in my simulation, ie. removing the coulombic interactions for a hydrogen atom.
Kind regards

Dries

Documentation clarification for u_klt matrix

Comments in the alchemical_analysis.py script and parser scripts state:

u_klt[k,m,t] is the reduced potential energy of snapshot t of state k evaluated at state m

But code has it pulling from the dUi=1,j columns in the .xvg files so it appears u_klt is a matrix of potential energies differences

mbar.getFreeEnergy differences returns an unexpected number of arguments

Hi,

The script failed with the following messages:

[ snip... ]
Final dimensionless free energies
f_k =
[ 0. 4.55294129 8.48088947 10.41122752 11.36350297
11.95779766 11.98970287 9.30249125 7.6417943 ]
MBAR initialization complete.
Traceback (most recent call last):
File "/data/user/sw/alchemical-analysis/alchemical_analysis.py", line 1130, in
Deltaf_ij, dDeltaf_ij = estimatewithMBAR(u_kln, N_k, P.relative_tolerance, regular_estimate=True)
File "/data/user/sw/alchemical-analysis/alchemical_analysis.py", line 252, in estimatewithMBAR
(Deltaf_ij, dDeltaf_ij, theta_ij) = MBAR.getFreeEnergyDifferences(uncertainty_method='svd-ew')
ValueError: need more than 2 values to unpack

Could you please provide a solution to this? I would appreciate it. Here is the command line used to invoke the command:
$ python /data/user/sw/alchemical-analysis/alchemical_analysis.py -o analysis -p Ethanol.dhdl -s 0 -t 300 -v

I can provide the data I am using as well (its small, only 50 ps worth of data); if you need it, where should I put it so you can access it? Thanks for your help. Here is the header for lambda = 0:

This file was created Sun May 24 08:22:44 2015

Created by:

GROMACS: mdrun_mpi_d, VERSION 5.0.5 (double precision)

Executable: /usr/local/gromacs/bin/mdrun_mpi_d

Library dir: /usr/local/gromacs/share/gromacs/top

Command line:

mdrun_mpi_d -multi 9 -replex 1000 -nex 100000 -deffnm Ethanol. -dhdl Ethanol.dhdl.

mdrun_mpi_d is part of G R O M A C S:

Groningen Machine for Chemical Simulation

@ title "dH/d\xl\f{} and \xD\f{}H"
@ xaxis label "Time (ps)"
@ yaxis label "dH/d\xl\f{} and \xD\f{}H (kJ/mol [\xl\f{}]\S-1\N)"
@type xy
@ subtitle "T = 300 (K) \xl\f{} state 0: (coul-lambda, vdw-lambda) = (0.0000, 0.0000)"
@ view 0.15, 0.15, 0.75, 0.85
@ legend on
@ legend box on
@ legend loctype view
@ legend 0.78, 0.8
@ legend length 2
@ s0 legend "Total Energy (kJ/mol)"
@ s1 legend "dH/d\xl\f{} coul-lambda = 0.0000"
@ s2 legend "dH/d\xl\f{} vdw-lambda = 0.0000"
@ s3 legend "\xD\f{}H \xl\f{} to (0.0000, 0.0000)"
@ s4 legend "\xD\f{}H \xl\f{} to (0.2000, 0.0000)"
@ s5 legend "\xD\f{}H \xl\f{} to (0.5000, 0.0000)"
@ s6 legend "\xD\f{}H \xl\f{} to (1.0000, 0.0000)"
@ s7 legend "\xD\f{}H \xl\f{} to (1.0000, 0.2000)"
@ s8 legend "\xD\f{}H \xl\f{} to (1.0000, 0.4000)"
@ s9 legend "\xD\f{}H \xl\f{} to (1.0000, 0.6000)"
@ s10 legend "\xD\f{}H \xl\f{} to (1.0000, 0.8000)"
@ s11 legend "\xD\f{}H \xl\f{} to (1.0000, 1.0000)"
@ s12 legend "pV (kJ/mol)"
0.0000 -28913.188 14.682067 8.8265472 0.0000000 2.9364134 7.3410336 14.682067 17.557088 21.936508 27.137445 32.781006 38.647377 1.6391496
0.2000 -28904.606 57.658838 -33.534042 0.0000000 11.531768 28.829419 57.658838 55.403482 57.872462 62.203406 67.342486 72.837012 1.6398901
0.4000 -28867.621 50.950735 -59.005125 0.0000000 10.190147 25.475367 50.950735 45.220515 46.121345 49.772451 54.682820 60.188709 1.6413282
0.6000 -28940.959 60.787208 -84.862704 0.0000000 12.157442 30.393604 60.787208 51.066626 49.699611 52.078965 56.290448 61.432133 1.6436758
0.8000 -29012.527 58.265203 -9.2633724 0.0000000 11.653041 29.132602 58.265203 58.815212 62.352019 67.351139 73.089061 79.191673 1.6479552
1.0000 -28913.114 52.140243 -56.167506 0.0000000 10.428049 26.070121 52.140243 46.579172 47.537230 51.377544 56.591160 62.473134 1.6529072
1.2000 -28995.072 54.985777 -60.363362 0.0000000 10.997155 27.492889 54.985777 48.927036 49.548740 52.984165 57.714439 63.059163 1.6569541
1.4000 -29047.081 53.638803 -28.937149 0.0000000 10.727761 26.819402 53.638803 51.889673 54.496837 58.859205 63.994079 69.471719 1.6605582
1.6000 -28980.959 48.438421 -49.952803 0.0000000 9.6876842 24.219211 48.438421 43.895239 45.418109 49.523529 54.831600 60.713221 1.6658269
1.8000 -29088.399 70.849208 -28.480365 0.0000000 14.169842 35.424604 70.849208 68.868205 71.363564 75.972778 81.635058 87.822310 1.6709045
2.0000 -29081.420 43.606014 -16.728997 0.0000000 8.7212029 21.803007 43.606014 43.122092 46.265671 51.235331 57.151142 63.548873 1.6758257
2.2000 -29100.045 62.089164 -16.022554 0.0000000 12.417833 31.044582 62.089164 62.068713 65.700084 70.937465 76.925393 83.257956 1.6783177
2.4000 -29091.456 69.023006 -59.653885 0.0000000 13.804601 34.511503 69.023006 63.206763 64.268652 68.302500 73.737728 79.850794 1.6788066
2.6000 -29137.602 79.838816 -55.844545 0.0000000 15.967763 39.919408 79.838816 74.157667 74.974791 78.801615 84.119636 90.190026 1.6785062
2.8000 -29137.617 114.95491 -76.349850 0.0000000 22.990981 57.477453 114.95491 106.58465 106.26599 109.59500 114.70361 120.70368 1.6781231

Sid

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.