paduagroup / fftool Goto Github PK

View Code? Open in Web Editor NEW

132.0 14.0 57.0 705 KB

Tool to build force field input files for molecular simulation

License: MIT License

Python 100.00%

molecular-dynamics force-field lammps packmol

fftool's Introduction

fftool

Agilio Padua

This is a Python tool to build force field input files for molecular dynamics simulations of systems composed of molecules, ions or extended materials.

fftool creates initial files for classical, fixed-charge molecular dynamics simulations. A force field database ionic liquids is available in CL&P. For polarizable force field and simulations, check the CL&Pol tools and database.

fftool: builds a simulation box and the corresponding force field for systems containing molecules, ions or extended materials. It requires the Packmol software to generate coordinates in the box. It outputs files in formats suitable for the LAMMPS, OpenMM, GROMACS or DL_POLY molecular dynamics packages.
tools/: utility scripts.
examples/: examples of molecule files and force field databases.

Requirements

Python
Packmol to pack molecules and materials in the simultion box.

Obtaining

Download the files or clone the repository:

git clone https://github.com/agiliopadua/fftool.git

Tutorial

These are instructions on how to build an force field files and an initial configuration for a system composed of molecules, ions or materials.

For each molecule, ion or fragment of a material prepare a file with atomic coordinates and eventually connectivity (covalent bonds). The formats accepted by this tool are .zmat, .xyz, .pdb or .mol, which are common formats in computational chemistry.

A .zmat file has the molecule name in the first line, followed by one empty line, then the z-matrix. See the examples directory and the Wikipedia entry Z-matrix(chemistry). Variables can be used in place of distances, angles and dihedrals. fftool infers connectivity (topology) from the z-matrix by default. In this case cyclic molecules require additional connect records to close rings. Improper dihedrals can be indicated by improper records. If a reconnect record is present, then connectivity will be guessed based on bond distances from the force field (see below). Below the z-matrix and the informations above, the name of a file with force field parameters can be supplied.

The XYZ file format .xyz contains atomic coordinates only. The name of a file with force field parameters can be given in the second line after the molecule name, and in this case connectivity is deduced from the bond lengths in the force field.

The PDB file format .pdb is widely used for proteins. The name of a file with force field parameters can be given on a COMPND record after the molecule name. fftool infers connectivity from the bond lengths in the force field (CONECT records are not read).

The MDL Molfile .mol file format contains a table with coordinates and also bonds. The name of a file with force field parameters can be given in the first line after the molecule name or in the third line. If the keyword reconnect is present after the force field filename, then connectivity will be deduced based on bond distances from the force field.

There are many tools (Open Babel, Avogadro, VESTA) to create file in the above formats. Manual editing of the files is usually necessary in order to match the atom names with those of the force field.
Use fftool to create an input file for packmol, which will use new _pack.xyz files with atomic coordinates for the components of your system. For help type fftool -h. For example, to build a simulation box with 40 ethanol and 300 water molecules and a density of 38.0 mol/L do:
```
 fftool 40 ethanol.zmat 300 spce.zmat -r 38.0
```
Alternatively, the side length of the the simulation box (here cubic) can be supplied in angstroms:
```
 fftool 40 ethanol.zmat 300 spce.zmat -b 20.0
```
Use packmol with the pack.inp file just created to generate the atomic coordinates in the simulation box:
```
 packmol < pack.inp
```
A difficult convergence may indicate that density is too high, so adjust density or box size if necessary. For more complex spatial arrangements of molecules and materials you can modify the pack.inp to suit your needs (see the Packmol documentation). Atomic coordinates for the full system are written to simbox.xyz.
Use fftool to create input files for LAMMPS (-l), OpenMM (-x), GROMACS (-g) ou DL_POLY (-d) containing the force field parameters and the coordinates of all the atoms (taken from simbox.xyz):
```
 fftool 40 ethanol.zmat 300 spce.zmat -r 38.0 -l
```
If no force field information was given explicitly in the molecule files, a default LJ potential with parameters zeroed will be assigned to atoms. No terms for bonds, angles or torsions will be created. This is suitable when working with non-additive, bond-order or other potentials often used for materials. The input files for MD simulations will have to be edited manually to include interaction potentials.

Deducing Bonds and Angles

When inferring connectivity from atomic coordinates, distances in the coordinates file are compared with equilibrium distances specified for bonds in the force field, and a tolerance of +/-0.25 angstrom is used to decide if a bond is created. So, the bond lengths in the conformation present in the molecule file must be close to those in the force field specification for those bonds to be included in the potential energy fonction of the system. The user is advised to check the number of bonds by creating a test system with the minimum of molecules.

Angles will be assigned to groups of three atoms i-j-k, with i-j and j-k bonded, if the value of the angle in the conformation from the molecule file is within +/-15 degrees of the equilibrium angle in the force field. If not, even if the atoms i-j-k are bonded, their angle will not be present in the final potential energy function, although topologically the angle is there. When running fftool to create a force field file (with -l, -x, -g or -d option) a warning message will show which such topological angles have been "removed" because they deviate too much from the equilibrium angles in the force field. This removal of angles avoids problems with atoms that have more than four ligands, such as S or P atoms with five or six ligands. Around these centers there are topological angles of 180 degrees to which no potential energy of bending is attributed in force fields. For example, in the octahedral PF6- anion there are two different values of F-P-F angles: twelve 90 degree angles between adjacent F atoms, and three 180 degree angles between opposite F atoms; only the twelve 90 degree angles contribute with a harmonic potential energy function in most force fields.

The tolerances for bond distances and angle values, 0.25 angstrom and 15 degrees, respectively, were chosen based on judgement. They can be set by editing the fftool source, namely the global variables BondTol and AngleTol. Use with care because spurious bonds and angles may be created if the tolerances are set too large.

Improper Dihedrals

Improper dihedrals are often used to increase the rigidity of planar atoms (sp2) and differ from proper dihedrals in how they are defined. A proper dihedral i-j-k-l is defined between bonded atoms i-j, j-k, and k-l and corresponds to torsion around bond j-k, the dihedral being the angle between planes i-j-k and j-k-l. An improper dihedral i-j-k-l is defined between bonded atoms i-k, j-k and k-l, therefore k is a central atom bonded to the other three. fftool assumes the central atom of the improper dihedral to be the third in the list. Often in force fields the same potential energy function is used both for proper and improper torsions.

If improper records are supplied in a molecule file (in .zmat format) then those improper dihedrals are read by fftool. Otherwise, the script will search for candidate improper dihedrals on all atoms with three bonds, whatever the input format. Warning messages will be printed if there are atoms with three bonds, and these messages can be ignored if the atoms in question are not centers of improper torsions. The user is advised to check the number and order of the atoms in the true improper dihedrals in the files created, by testing with a minimal system.

Periodic Boundary Conditions

In molecular systems the initial configuration will generaly not contain molecules crossing boundaries of the simulation box. A buffer distance of 1.5 angstrom is reserved at the box boundaries to avoid overlap of molecules from periodic images in the initial configuration, as explained in the packmol documentation (this empty space is added by fftool only for orthogonal boxes). So the user should be aware of this empty volume when choosing the size of the box.

For simulations with extended materials it is possible to create chemical bonds across boundaries. Option -p allows specification of periodic conditions along x, y, z or combinations thereof. It is important in this case to supply box dimensions using the option -b <l> for a cubic box, -b <lx,ly,lz> for a general orthogonal box, or -b <a,b,c,alpha,beta,gamma> for a general parallelepiped (triclinic box). An energy minimization step prior to the start of the MD simulation is highly recommended because fftool will leave no extra space near the boundaries and certain molecules may overlap with those of neighboring images.

The coordinates of the atoms of the material have to be supplied in .xyz format and prepared carefully so that distances across periodic boundaries are within the tolerance to identify bonds. The user is advised to check the number of bonds in the output files created.

It is important that only the material for which bonds are to be established across boudaries is supplied in .xyz format. The initial files for other molecules in the system should be in .zmat or .mol formats, which contain connectivity information. This is to avoid spurious bonds between atoms of the molecular species that may happen to be positioned too close to boundaries.

The pack.inp file will likely need manual editing in order to position the atoms of the material precisely.

Force Field File Format

The fftool script reads a database of molecular force field parameters in xml format (similar to the format used by OpenMM), or in the original .ff format described below. See the examples directory.

The ff2xml script converts from the original to the xml format.

.ff format

Blank lines and lines starting with # are ignored.

There are five sections, with headings ATOMS, BONDS, ANGLES, DIHEDRALS and IMPROPER. Under each section heading, registers concerning the different types of term in the force field are given.

ATOMS records describe, for each type of atom:

the non-bonded atom type used for intermolecular interactions (these types may differ in the charges or intermolecular potential parameters)
the bonded atom type used in intermolecular interactions (these types determine the intramolecular terms such as bonds, angles dihedrals)
the mass in atomic units
the electrostatic charge in elementary units
the non-bonded potential type, e.g. lj
potential parameters, namely Lennard-Jones sigma (angstrom) and epsilon (kJ mol-1)
```
  C3H   CT  12.011  -0.18   lj    3.50   0.27614
```

BONDS records describe covalent bonds between intramolecular atom types:

two bonded atom types
type of bond potential, e.g. harm for harmonic potential or cons for a constrained bond.
bond potential parameters, namely equilibrium distance (angstrom) and force constant in the form k/2 (x - x0)^2 (kJ mol-1 A-2)
```
  CT  CT   harm   1.529   2242.6
```

ANGLES records describe valence angles between intramolecular atom types:

three bonded atom types, in which the central atom is bonded to the other two, e.g. i-j and j-k are bonded.
type of angle potential, e.g. harm for harmonic potential or cons for a constrained angle.
angle potential parameters, namely equilibrium angle (degrees) and force constant in the form k/2 (x - x0)^2 (kJ mol-1 rad-2)
```
  HC  CT  CT   harm   110.7   313.8
```

DIHEDRALS records describe torsion angles between intramolecular atom types:

four bonded atom types, in which atoms i-j, j-k, k-l are bonded.
type of dihedral potential, e.g. opls for OPLS cosine series with four terms.
dihedral potential parameters, with the coefficients in the form V_n/2 (1 +/- cos(n phi)) (kJ mol-1).
```
  CT  CT  CT  CT   opls    5.4392   -0.2092    0.8368    0.0000
```

IMPROPER records describe improper dihedral angles between intramolecular atom types:

four bonded atom types, in which atoms i-k, j-k, k-l are bonded.
type of dihedral potential, e.g. opls for OPLS cosine series with four terms.
dihedral potential parameters, with the coefficients in the form V_n/2 (1 +/- cos(n phi)) (kJ mol-1).
```
  CA  CA  CA  HA   opls    0.0000    9.2048    0.0000    0.0000
```

References

Packmol: L. Martinez et al. J Comp Chem 30 (2009) 2157, DOI: 10.1002/jcc.21224
LAMMPS: S. Plimton, J Comp Phys 117 (1995) 1, DOI: 10.1006/jcph.1995.1039
OpenMM: P. Eastman, J. Swails, J. D. Chodera, R. T. McGibbon, Y. Zhao, K. A. Beauchamp, L.-P. Wang, A. C. Simmonett, M. P. Harrigan, C. D. Stern, R. P. Wiewiora, B. R. Brooks, and V. S. Pande. PLOS Comp. Biol. 13 (2017) e1005659, DOI: 10.1371/journal.pcbi.1005659
GROMACS: H.J.C. Berendsen, D. van der Spoel, R. van Drunen, Comp Phys Commun, 91 (1995) 43, DOI: 10.1016/0010-4655(95)00042-E
DL_POLY: I.T. Todorov and W. Smith, Daresbury Lab.

fftool's People

Contributors

Stargazers

Watchers

Forkers

rocli izxle iurisegtovich zhujieli shrubag trollchu turneradam like201707 zhongxiang117 plin1112 ananthgr tjkartha georgeathust kateryna-goloviznina longkunxuluke tommason14 warmshawn lorenaav yongjiguan dezhao-huang yhuang-thu rberthin naterriver peachrl chengjixian sheng-bi mrat1618 shaocheng16 hrsafavi eastsheng kmmukut naokicc lorenzofollador whitehorseprince abdulelahalshehri christianzhao kawai-ryotaro rnaimehaom jimoo068 zhangdft kylincaster xiaoyu-wang-stone lonxunquantum park-hyungshick similt xcnmvsdj mondesinxi tanshendong2023 chrinide mcancade nancysnals alinacao2000 badatoco zhouwinnie askforarun currymultisci

fftool's Issues

LAMMPS input order -> tail correction lost after restart

Dear Agilio Padua,

I have been using fftool for preparing my MD simulations with LAMMPS. I recently fell prey to an unexpected effect due to the order of directives in the LAMMPS input scripts generated by fftool.
If one uses pair_style hybrid, the pair_modify settings (e.g., tail correction) don't seem to be stored in the binary restart file. This will lead to an unexpected, yet significant, drop in the density in a simulation.

I have submitted the issue to the LAMMPS-users mailing list (see link below), but it is yet unclear whether this is going to be resolved internally in LAMMPS. To avoid any possible issues and confusion I would suggest changing the order of the restart and pair_style/pair_modify directives in the LAMMPS input generation in fftool.

If you want I can create a corresponding pull request.

Cheers,

Frank

inquiry to the LAMMPS-users mailing list: https://sourceforge.net/p/lammps/mailman/message/35456160/

pdb files for systems with more than 99999 atoms

Hello all,

I have generated input files for OpenMM using fftool and my system contains more than 100000 atoms. When running OpenMM I get an error "ValueError: Misaligned residue name" and it points to the line of atom 100000.
The error is solved by manually removing one space from the front of the residue name for all residues after 99999.
Is it possible to include this on the code?

Thank you!

Units conversion: GROMACS to LAMMPS

Good morning,

I am having some issue in understanding the units conversion from the GROMACS .itp file into the LAMMPS .ff one. Please correct me if I am wrong.
GROMACS adopts the following units:

sigma [nm] epsilon [kJ/mol] b_0 [nm] k_b [kJ mol-1 nm-2]

whereas LAMMPS, when units real key is used, has the following units:

sigma [Å] epsilon [kcal/mol] b_0 [Å] k_b [kcal mol-1 Å-2]

To give you an example the .itp file is like:

[ atomtypes ] opls_800 C800 12.0110 0.000 A 3.50000E-01 2.76144E-01, where [ atoms ] 1 opls_800 1 UNK C00 1 -0.0487 12.0110

Whereas .ff looks like:

ATOMS C00 C00 12.011 -0.0487 lj 3.50 0.27614

Here I see the conversion from nm to Å, but I don't see the conversion from kJ/mol to kcal/mol. I am quite sure that here I am missing something, but I don't what I am missing.

Furthermore the equation used by GROMACS to define the bonded interaction uses a 1/2 term as prefactor, whereas the one used by LAMMPS, when the bond_style harmonic key is used, don't. Which means that all the k_b terms should be divided by 200 (if we don't consider the kJ to kcal conversion) when converting from GROMACS ([kJ mol-1 nm-2]) to LAMMPS real units ([kcal mol-1 Å-2]). Whereas I get:

GROMACS .itp
[ bonds ] 2 1 1 0.1410 267776.000

LAMMPS .ff
BONDS O01 C00 harm 1.410 2677.8

I am probably missing something since the simulations I run give reasonable results, but I just wanted to be sure about why especially since I need to write a customized .ff file by hand and I don't want to use erroneous values/units.

Thanks.

Matteo

output to Amber

Is there a way for fftool to output to Amber? I tried to output to GROMACS first and then used Parmed to convert GROMACS topology to Amber but got the error that it couldn't identify the file format.

Parameters for cyclic carbonates

Hi, I was trying to simulate some carbonates, specifically propylene carbonate and ethylene carbonate and other not cyclic, but I can't find the parameters in the force field so that the molecules are neutral. The cyclic residuals with the respective parameters used among those present in the force field are:

  O_2                 O_2
  ||                  ||
  C_2               C_2
 /   \               /   \
OS   OS     OS   OS
 \   /               \   /
 CT-CT        CT-CT
         \           
          CT

propylencarbonate ethylencarbonate

Thank you for your time.

Unable to read coordinates files of custom molecule

I have tried to read a molecule with fftool and in several data formats but I encounter several issues:
The molecule is TFSI:

16

C          6.03034        4.09220        5.48886
F          5.37202        5.30065        5.36546
F          6.39135        3.89541        6.80769
F          7.16406        4.10454        4.69920
N          3.42469        2.54403        5.86881
C          0.86033        0.89196        5.97521
F          0.03319        0.14152        5.16196
F          0.14535        1.94693        6.50866
F          1.34449        0.09819        6.99698
S          4.93349        2.75026        4.95267
S          2.24665        1.53550        4.99964
O          5.62532        1.37952        5.01298
O          4.53998        2.86277        3.47122
O          1.78341        2.37997        3.80211
O          3.05533        0.41678        4.32326
H          3.00037        3.48797        6.01118

I created a .pdb and a .gzmat extension with obabel and a .zmat extension with molden. All give different errors.

$ fftool 10 files/TFSI.xyz -r 10
> density 10.000 mol/L  volume 1660.6 A^3
> molecule_file      species           nmol force_field      nbond source  charge
> Traceback (most recent call last):
  File "/Users/marcodigennaro/WORK/external_packages/fftool/fftool", line 2816, in <module>
    main()
  File "/Users/marcodigennaro/WORK/external_packages/fftool/fftool", line 2768, in main
    spec.append(mol(zfile, connect, box))
  File "/Users/marcodigennaro/WORK/external_packages/fftool/fftool", line 615, in __init__
    self.fromxyz(filename, connect, box)
  File "/Users/marcodigennaro/WORK/external_packages/fftool/fftool", line 882, in fromxyz
    self.name = tok[0]            # molecule name
IndexError: list index out of range

$ fftool 10 files/TFSI.pdb -r 10
> density 10.000 mol/L  volume 1660.6 A^3
> molecule_file      species           nmol force_field      nbond source  charge
>  warning: no force field for TFSI
> Traceback (most recent call last):
  File "/Users/marcodigennaro/WORK/external_packages/fftool/fftool", line 2816, in <module>
    main()
  File "/Users/marcodigennaro/WORK/external_packages/fftool/fftool", line 2772, in main
    .format(zfile, spec[i].name, spec[i].nmol, spec[i].ffile,
TypeError: unsupported format string passed to NoneType.__format__

$ fftool 10 files/TFSI.zmat -r 10
> density 10.000 mol/L  volume 1660.6 A^3
> molecule_file      species           nmol force_field      nbond source  charge
> Traceback (most recent call last):
  File "/Users/marcodigennaro/WORK/external_packages/fftool/fftool", line 2816, in <module>
    main()
  File "/Users/marcodigennaro/WORK/external_packages/fftool/fftool", line 2768, in main
    spec.append(mol(zfile, connect, box))
  File "/Users/marcodigennaro/WORK/external_packages/fftool/fftool", line 611, in __init__
    self.fromzmat(filename, connect)
  File "/Users/marcodigennaro/WORK/external_packages/fftool/fftool", line 696, in fromzmat
    self.zmat2cart(z)
  File "/Users/marcodigennaro/WORK/external_packages/fftool/fftool", line 754, in zmat2cart
    theta = math.acos(delx / math.sqrt(delx*delx + dely*dely))
ZeroDivisionError: float division by zero

How to get the right zmat format? In the example folder some files (eg. ethane) have leading atom index:

ethane

1  C3H
2  C3H  1  rCC
3  H    2  rCH  1  aHCC
4  H    1  rCH  2  aHCC  3   d1
5  H    2  rCH  1  aHCC  4   d2
6  H    1  rCH  2  aHCC  5   d1
7  H    2  rCH  1  aHCC  6   d2
8  H    1  rCH  2  aHCC  7   d1
...

others (eg. benzene) do not:

benzene

CA
CA  1  rCC
CA  2  rCC  1  120.00
CA  3  rCC  2  120.00  1   0
CA  4  rCC  3  120.00  2   0
CA  5  rCC  4  120.00  3   0
HA  1  rCH  2  120.00  3 180
HA  2  rCH  3  120.00  4 180
HA  3  rCH  4  120.00  5 180
HA  4  rCH  5  120.00  6 180
HA  5  rCH  6  120.00  1 180
HA  6  rCH  1  120.00  2 180
...

What is the correct format and how to obtain it?
In my case, I don't know the atom types in advance (TFSI is just and example).
Thank you

fftool: error: the following arguments are required: infiles

Thanks, Now, Solved!

Newest commit creates incorrect pdb files with ilff

Master created config.pdb on the left, fftool-02/22/2019 on the right (runs in GROMACS directly)

PDB warning message

Do we need a warning message when generating an input file for LAMMPS and not for OpenMM?

  warning: atom name SiO1a too long for pdb in siosicmim+
  warning: atom name HCSi1 too long for pdb in siosicmim+
...

xyz and ff file: improper keyword isnt detected

I am using your fftool to convert a xyz file to the lammps input files and ran into some trouble with the IMPROPER keyword in the ff file, it seems the code doesnt read it. Therefor the impropers get added to the list of dihedrals. Also the potential keyword isnt read for dihedral and defaults to opls. I have added the xyz and ff files for reference.

xyz.txt
ff.txt