Skip to content

cecimarques/gff_lammps

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

25 Commits
 
 
 
 
 
 
 
 
 
 
 
 

Repository files navigation

Hello :)

Here you can find a python code (namely code.py) built for finding the OPLS-AA force field parameters for any given compound in GROMACS .itp files and outputting them in the LAMMPS format considering "real" units. The code can be adapted for any force field, with OPLS-AA being merely an example, AND was designed to work in conjuncture with the 1-molecule LAMMPS data file that is generated by LigParGen (https://zarbi.chem.yale.edu/ligpargen)* from the compound's SMILES. Sample files specifically for pentanol are given in the "Sample-input" directory. These go together with code.py as it is given here.
* LigParGen does output the OPLS-AA potential parameters in the LAMMPS data file, but (in my experience!!!) these were found to be bogus in the case of many compounds compared to what I would expect from the original publication for example *and* the charges sometimes dont sum to 0 for compounds that should be charge neutral.

The files ffnonbonded.itp, ffbonded.itp and atomtypes.itp found in this directory together with the python code were directly taken from the official GROMACS distribution (version 2024.4, date released 31/10/2024) github page for the OPLS-AA force field (https://github.com/gromacs/gromacs/tree/main/share/top/oplsaa.ff). The two former were used to build input files required by the python code while the latter is necessary as "guidance" to assign suitable "OPLS-AA atom types" within the nomenclature used by GROMACS (see discussion below - the code requires *you* to choose duly the "OPLS-AA atom types" for your compound).

INPUT FILES REQUIRED BY THE PYTHON CODE (8 in total):
- ffnonbonded : consists of a straightforward copy-paste of the [ atomtypes ] section of the OPLS-AA ffnonbonded.itp file.
- ffbonds : consists of a straightforward copy-paste of the [ bondtypes ] section of the OPLS-AA ffbonded.itp file.
- ffangles : consists of a straightforward copy-paste of the [ angletypes ] section of the OPLS-AA ffbonded.itp file.
- ffdihedrals : consists of a straightforward copy-paste of the [ dihedraltypes ] section of the OPLS-AA ffbonded.itp file.
- bonds: file containing the content found in the "Bonds" section of the 1-molecule LAMMPS data file output by LigParGen (edited so that the columns are separated by "\t").
- angles: file containing the content found in the "Angles" section of the 1-molecule LAMMPS data file output by LigParGen (edited so that the columns are separated by "\t").
- dihedrals: file containing the content found in the "Dihedrals" section of the 1-molecule LAMMPS data file output by LigParGen (edited so that the columns are separated by "\t").
- impropers: file containing the content found in the "Impropers" section of the 1-molecule LAMMPS data file output by LigParGen (edited so that the columns are separated by "\t").

PS: Note that sometimes bonded potential parameters information concerning specific directives are found split throughout the ffbonded.itp file. In the latter case, you may choose to ignore information that is not relevant for the kind of compound you are going to use this code for (e.g. having lines for dihedrals that specifically concern some amino acids is not relevant if using the code to get potential parameters for organic solvents).

Comments are presented along the code to assist the reader in understanding it, and an overview of the code can be found below:

 -------------------------------------------------------- CODE OVERVIEW ----------------------------------------------------
a) The first step is to assign the "OPLS-AA atom types" within GROMACS nomenclature for the atoms in a molecule of the compound of interest using the "atomtypes.atp" file. This should be doable under the assistance of the comment lines found in that file. These "OPLS-AA atom types" for each atom are then declared in the python code manually in ascending order with respect to the LAMMPS atom-IDs (given in the 1-molecule LigParGen LAMMPS data file). The user may use a visualization software (e.g. OVITO) to better grasp the correspondence between the LAMMPS atom-IDs and where the atom is within the molecule.
b) The code is going to read the "ffnonbonded" file and search for the lines that concern the "OPLS-AA atom types" you input and save their values of "charge", "epsilon", "sigma" and their "atom type within OPLS-AA nomenclature" (i.e. name that appears in the second column of the ffnonbonded.itp file). The epsilon and sigma values are going to be output in the kernel in a format that can be directly copied/pasted into the "Pair Coeffs" section of the LAMMPS data file. There will be one (and only one, since no repetition of info for a given "OPLS-AA atom type" should exist in the ffnonbonded.itp file by GROMACS!) line printed for each LAMMPS atom ID. The possibility of copying/pasting directly is unlocked by the fact that LigParGen outputs a 1-molecule LAMMPS data file that has a distinct "LAMMPS atom type" for each atom in the molecule. An array containing the charges is going to be saved into a file and can be easily plugged into 3rd column of the "Atoms" section (atom_style format full) of the LAMMPS data file. There is also a built-in "safety check" that sums the charges of the atoms composing a molecule (should be 0 if the molecule is neutral): thet net charge of the molecule will be printed and you may check it in the kernel if you want.
c) The code reads the "bonds", "ffbonds", "angles", "ffangles", "dihedrals", "ffdihedrals" files and uses them to find the potential parameters underlying each bond, angle and dihedral declared in the sections of the 1-molecule LAMMPS data file gotten from LigParGen. This is done by associating the LAMMPS atom IDs for each atom with the "atom type within OPLS-AA nomenclature" and later finding the line in the "ffbonds", "ffangles" and "ffdihedrals" files that contain the parameters for the potential concerning the specific sequence of "atom types within OPLS-AA nomenclature". This information is also output in the kernel in a format that allows direct copying/pasting in the "Bond Coeffs", "Angle Coeffs" and "Dihedral Coeffs" LAMMPS data file. Note that, similarly to the case of the "epsilon" and "sigma" LJ parameters in the "Pair Coeffs" section, this copying/pasting is only possible thanks to the fact that LigParGen assigns a unique LAMMPS bond/angle/dihedral type to each LAMMPS bond/angle/dihedral ID. IMPORTANT!: the copying/pasting however needs to be done with caution because (1) the code outputs twice the parameters found in the lines of the ffbond, ffangle and ffdihedral that contain bonds, angles and dihedrals that are formed by a palindrome sequence of "atom type within OPLS-AA nomenclature" and (2) the code does not check for bonds, angles and dihedrals formed by a wildcard atom type, which then need to be searched manually in the .itp file. Thus, the copying/pasting needs to be done mindfully (meaning deleting repeated lines AND checking if indeed parameters for all bonds, angles and dihedrals were found). 

The .py code takes care of all necessary unit conversions required in the process of transforming GROMACS potential parameters to LAMMPS potential parameters for a LAMMPS input script that contemplates units "real". It also accounts for the differences in the potential form as implemented in GROMACS and LAMMPS (e.g., sometimes a pre-factor of 1/2 exists). More specifically, the parameters are set to go along with an input script that cherishes the following setup:
pair_style      lj/cut/coul/long 11 11
pair_modify     mix geometric tail yes
bond_style      harmonic
angle_style     harmonic
dihedral_style  opls
improper_style  cvff
Further specific setup:
special_bonds   lj 0 0 0.5 coul 0 0 0.5
kspace_style    pppm 1e-6

d) Finally, for impropers, the code simply outputs the list of "atom type within OPLS-AA nomenclature" underlying each improper declared below the "Impropers" section in the 1-molecule LAMMPS data file gotten from LigParGen. It is responsibility of the user to see, based on the sequence of atom types underlying each improper, if it is necessary to set an improper potential for any of the impropers. Furthermore, it is worth mentioning that GROMACS advises implementing improper potentials as dihedral potentials and the corresponding potential form prescribed can be found in the original ffbonded.itp file in the GROMACS tree. You may choose to set them as improper potentials instead within LAMMPS (up to you)!
 -----------------------------------------------------------------------------------------------------------------------------

Hope you can find it useful !

- Cecilia Alvares and Gabriele Sosso

About

Python code to get force field parameters from GROMACS .itp files

Resources

Stars

Watchers

Forks

Releases

No releases published

Packages

No packages published

Languages