GMXPBSA 2.1 is a user-friendly suite of Bash/Perl scripts for streamlining MM/PBSA calculations on structural ensembles derived from GROMACS trajectories, to automatically calculate binding free energies for protein–protein or ligand–protein complexes [R.T. Bradshaw et al., Protein Eng. Des. Sel. 24 (2011) 197–207]. GMXPBSA 2.1 is flexible and can easily be customized to specific needs and it is an improvement of the previous GMXPBSA 2.0 [C. Paissoni et al., Comput. Phys. Commun. (2014), 185, 2920–2929]. Additionally, it performs computational alanine scanning (CAS) to study the effects of ligand and/or receptor alanine mutations on the free energy of binding. Calculations require only for protein–protein or protein–ligand MD simulations. GMXPBSA 2.1 performs different comparative analyses, including a posteriori generation of alanine mutants of the wild-type complex, calculation of the binding free energy values of the mutant complexes and comparison of the results with the wild-type system. Moreover, it compares the binding free energy of different complex trajectories, allowing the study of the effects of non-alanine mutations, post-translational modifications or unnatural amino acids on the binding free energy of the system under investigation. Finally, it can calculate and rank relative affinity to the same receptor utilizing MD simulations of proteins in complex with different ligands. In order to dissect the different MM/PBSA energy contributions, including molecular mechanic (MM), electrostatic contribution to solvation (PB) and nonpolar contribution to solvation (SA), the tool combines two freely available programs: the MD simulations software GROMACS [S. Pronk et al., Bioinformatics 29 (2013) 845–854] and the Poisson–Boltzmann equation solver APBS [N.A. Baker et al., Proc. Natl. Acad. Sci. U.S.A 98 (2001) 10037–10041]. All the calculations can be performed in single or distributed automatic fashion on a cluster facility in order to increase the calculation by dividing frames across the available processors. This new version with respect to our previously published GMXPBSA 2.0 fixes some problem and allows additional kind of calculations, such as CAS on single protein in order to individuate the hot-spots, more custom options to perform APBS calculations, improvements of speed calculation of APBS (precF set to 0), possibility to work with multichain systems (see Summary of revisions for more details). The program is freely available under the GPL license. New version program summaryProgram title: GMXPBSA 2.1Catalogue identifier: AETQ_v1_1Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AETQ_v1_1.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: GNU General Public License, version 3No. of lines in distributed program, including test data, etc.: 188453No. of bytes in distributed program, including test data, etc.: 5381839Distribution format: tar.gzProgramming language: Bash, Perl.Computer: Any computer.Operating system: Linux, Unix OS.RAM: 2147483648 bytesSupplementary material: A table of added keywords is available.Classification: 3.Catalogue identifier of previous version: AETQ_v1_0Journal reference of previous version: Comput. Phys. Comm. 185 (2014) 3016External routines: Any APBS version (http://www.poissonboltzmann.org/apbs/) and GROMACS version 4.5 or later installations (http://www.gromacs.org). Optionally LaTeX.Does the new version supersede the previous version?: YesNature of problem:The Molecular Mechanics (MM) data (Lennard-Jones and Coulomb terms) and the solvation energy terms (polar and nonpolar terms respectively) from an ensemble of structures derived from GROMACS molecular dynamics simulation trajectory are calculated. These calculations are performed for each single component of the simulated complex, including protein and ligand. In order to cancel out artifacts an identical grid setup for each component, including complex, protein and ligand, is required. Statistical analysis on the extracted data and comparison with wild-type complex in the case of either computational alanine scanning or calculations on a set of simulations is performed. Possible outliers in the frames extracted from the simulations during the binding free energy calculations are evaluated.Solution method:The tool combines the freely available programs GROMACS and APBS to: 1.extract frames from a single or multiple complex molecular dynamics (MD) simulation, allowing comparison between multiple trajectories;2.split the complex frames into the single components including complex, protein and ligand;3.calculate the Lennard-Jones and Coulomb energy values (MM terms);4.calculate the polar solvation energy values using the implicit solvation Poisson–Boltzmann model (PB);5.calculate the nonpolar solvation energy value based on the solvent accessible surface area (SASA);6.combine all the calculated terms into the final binding-free energy value;7.repeat the same procedure from point 1 to 6 for each simulation in the case of computational alanine scanning (CAS) or simultaneous comparison of different Mds.Reasons for new version:We had a lot of feedback from our previous GMXPBSA 2.0 publication suggesting improvements. There were also some bugs that required to be fixed urgently in order to keep the program working properly.Summary of revisions:In the update version of GMXPBSA 2.1 we include the following changes: 1.Fixed bug related to the use of amino acids with id>1000.2.Fixed bug related to the computing of the percentage of failed APBS-jobs.3.Added the string ‘export LC_NUMERIC = “en_us.UTF-8”’ in each of the three gmxpbsa*.sh scripts to avoid errors related to locale environment variables.4.Fixed a bug related to the CAS mutation in the cases in which: 1.- no chains information is available from xtc/tpr2.- the ligand is positioned before the receptor in the pdb5.Added final message with proper references.6.Added the possibility of checking the GMXPBSA version with the command “sh gmxpbsa[0,1,2].sh -h”.7.Included the possibility to use “multichain”. (Added keyword “multichain”, set by default to “n”). Useful if in the pdbs extracted from the trajectory there is more than one chain. The option “multichain” must be used ONLY if the string TER is not present at the end of each chain in the comp/receptor pdb files.8.Added the possibility of performing ΔG CAS calculations on a single protein. Added keyword “protein_alone”, set by default to “n”, and keyword “itp_protein” that must be set only in the case of “use_topology = y”. When using the option “protein_alone” in the xtc file provided to GMXPBSA the periodic boundary condition must be treated so as to make the broken molecules whole (i.e. trjconv -pbc whole).9.Modified default values for keywords: coulomb from “coul” to “gmx”, precF from “1” to “0”, linearized from “n” to “y”.10.Added the possibility of using custom name of trajectories and binary file, xtc and tpr respectively format.Many keywords have been added to allow more flexibility with the APBS options in the input file and to set the option for the cluster. See Table 1 of supplementary material for a list of all the added keywords.Restrictions:Input format files compatible with GROMACS engine 4.5 and later versions. Availability of the force field or of the topology files.Running time:On a single core Lennard-Jones, Coulomb and nonpolar solvation term calculations require a few minutes. The time required for polar solvation term calculations depends on the system size.