Generalized Amber Force Field (GAFF) for polymer electrolyte systems
This is a tutorial to build initial configurations and force field parameters required for the Classical Molecular dynamics simulation using GROMACS software of polymer electrolyte system PEO-LiTFSI. This tutorial is based on the jupyter notebook, written using different open-source tools. Before building the PE system, one has to answer the following questions:
| Questions | In this tutorial |
|---|---|
| What is the monomer? | -COCOCOCOC- |
| No. of monomers in the polymer chain? | 20 |
| What is the end group? | -CH_3 |
| What is the salt? | LiTFSI |
| No. of polymers and salt ions | 100,20 |
This tutorial assumes the user is familiar with Python (jupyter notebook) and MD simulations using GROMACS. The prerequisites for running the jupyter notebook are:
- GROMACS
- ambertools
- antechamber
- prepgen
- tleap
- Python packages
- jupyter
- ASE
- nglview
- intermol
Run the initial cell to import python packages required for the tutorial and create a working directory.
import numpy as np
import os
import nglview as ngl
from ase.io import read, write
WDIR='/home/harish/GroPolBul/tutorial/PEO-LiTFSI'
%mkdir {WDIR}
%cd {WDIR}The monomer or short polymer chain can be read from any structural file formats such as *.pdb, *.xyz, and *.gro and then "antechamber" can be used to parameterize the monomer.
#Use short polymer with repeat group, head and tail group
#repeat group: COCCOCCOCCOCCOC
#head and tail group: C
#Use antechamber (ambertools) to optimize and prameterise with AM1-BCC charges
#antechamber options help:
#Usage: antechamber -i input file name
# -fi input file format
# -o output file name
# -fo output file format
# -c charge method
# -nc net molecular charge (int)
# -rn residue name
!antechamber -i PEO_initial.pdb -fi pdb -o PEO.ac -fo ac -at gaff -an y -c bcc -nc 0 -rn PEO The monomer can also be visualized using nglview, and atom indexes are labeled, which is helpful to define the HEAD,CHAIN, and TAIL parts of the monomer.
#Visualize the molecule
mol=read('PEO_initial.pdb')
vi=ngl.show_ase(mol);vi.add_label(radius=2,color='black',label_type='atomindex')
vi
{width=100px}
This is the most complicated step of this tutorial as it required to choose the two different atom indexes:
- Head/Tail atom index numbers
- The atom indexes needed to be omitted near Head/Tail ids.
#Defining CHAIN HEAD and TAIL in the monomer directory
ac=open('PEO.ac',mode='r') #Reading .ac file
[next(ac) for _ in range(2)] #Skipping first two lines of text
l=ac.readlines() #Reading lines
#Breaking the parts of monomer to CHAIN, HEAD and TAIL
#Atom index where head and tail of monomer; Check from above ngl view of mol
#Change these values accordingly
head_id=1
tail_id=7
head_omit=[3, 23, 24, 25] #Atoms to omit near head
tail_omit=[8, 32, 33, 34] #Atoms to omit near tailUsing "prepgen" we will generate three *.prepi files
- PEO.prepi : CHAIN
- HPT.prepi : HEAD part
- TPT.prepi : TAIL part
Using these files, we can build the polymer chain.
Here we need to set the variable "n_mono_repeat" which saves number of monomer units in the short polymer chain, and need to choose the variable "n_mono_pol" i.e. number of monomers in the desired polymer. Next "tleap" is used to prepare the desired single polymer chain.
#Build the desired homopolymer using above prepi files
%mkdir polymer
n_mono_repeat=5
n_mono_pol=20
repeat=" ".join(['PEO'] * int(int(n_mono_pol/n_mono_repeat)-2))
print('HPT '+str(repeat)+' TPT')
#Write the tleap input file to combine the preparatory files and build polymer chain
tleap=open('polymer/PEO_tleap.in','w+')
tleap.write('''source leaprc.gaff
loadamberprep PEO.prepi
loadamberprep HPT.prepi
loadamberprep TPT.prepi
mol = sequence {HPT '''+str(repeat)+''' TPT}
savepdb mol polymer/PEO_'''+str(n_mono_pol)+'''mer.pdb
saveamberparm mol polymer/PEO_'''+str(n_mono_pol)+'''mer.prmtop polymer/PEO_'''+str(n_mono_pol)+'''mer.inpcrd
quit''')
tleap.close()
!tleap -s -f polymer/PEO_tleap.in > polymer/PEO_tleap.outThis will generate *.prmtop and *.inpcrd files, these are coordinates and parameter files for AMBER. These files can be converted from AMBER to GROMACS using "intermol".
#Converting AMBER to GROMACS using "intermol".
!python /home/harish/InterMol/intermol/convert.py --amb_in polymer/PEO_{n_mono_pol}mer.inpcrd polymer/PEO_{n_mono_pol}mer.prmtop --gromacsFinally visualize the polymer using nglview.
{width=100px}
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent