Modeling Oxidized Cellulose Crystals
Published:
Finally figured out the modeling of oxidized cellulose crystals with the CHARMM force field. The approach is similar for other force fields. LAMMPS’s built-in modeling capabilities are really powerful. Mark this down — it deserves a thorough write-up!
Overview
The most common oxidized form of cellulose crystals is TEMPO-oxidized cellulose. Specifically, the hydroxymethyl groups on the cellulose crystal surface are oxidized to carboxyl groups, with an oxidation degree of 0.1~0.5 (note: the degree of sulfonation is generally 0.1~0.2). For the CHARMM force field, some researchers have already achieved cellulose crystal modeling (cellulose-toolkit), but due to the complexity of model construction, studies on surface oxidation are relatively limited. (Existing research has used the OPLS-AA force field: Biomacromolecules 2020, 21, 3069−3080)
Basic Approach
Generally, force field parameters are included inside the data file for polymer systems.
Two reasons: 1. Polymer systems have many parameters, making it inconvenient to write them all in the input file. 2. The pair_coeff between different atom types is calculated according to the mixing rule (pair_modify), so pair_coeff values differ under different mixing rules. It is simpler to directly write the pair_coeff for individual atom types in the data file and declare the mixing rule in the input file. Alternatively, force field parameters can be included directly using the include command, allowing force field parameters and topology to be read separately. Two notes: 1. Although force field parameters and topology are read separately, the topology must be consistent with the corresponding force field — for example, you cannot read CHARMM force field parameters while using an OPLS-AA topology, because the atom types differ between force fields. 2. Force field parameters can have redundancy! This is important — having extra unused parameters is not a problem. With these two points in mind, modeling oxidized cellulose crystals essentially consists of three parts:
- 1 Obtain the Coeff parameter file for oxidized cellulose crystals
- 2 Modify the topological structure of the cellulose crystal
The first step can be done by adding force field parameters to the cellulose crystal data file, then using LAMMPS’s built-in write_coeff command to obtain the desired parameter file. (Making good use of LAMMPS’s own modeling capabilities can solve many problems.) Adding force field parameters can be done via a Python script. Note that the input data file must be a standard LAMMPS output data file.
# -*- codeing = utf-8 -*-
# @Time :2022/9/9 23:34
# @Author : 得意喵 ~
# @File : addcoeff.py.py
# @Software: PyCharm
import pandas as pd
import numpy as np
import re
#define functions
def findcoeff(data,coeff):
for i in range(len(data)):
if re.match(coeff,data[i]):
return i
def coeffformat(n,coeff0):
return str(n)+' '+coeff0+'\n'
def getcoeffindex(coeff0):
s=coeff0.split()
return int(s[0])
def addcoeffs(data,n,coeff,n0):
for i in range(len(coeff)):
data.insert(n+i,coeffformat(n0+i+1,coeff[i]))
return data
#open files
f=open('try.data','r')
data=f.readlines()
if data[10]=='\n':
data.insert(10,'0 impropers\n')
data.insert(11,'0 improper types\n')
data.insert(findcoeff(data,'Atom')-1,'\n')
data.insert(findcoeff(data,'Atom')-1,'Improper Coeffs\n\n')
#types and coeffs
new_atomtypes=2
new_bondtypes=2
new_angletypes=5
new_dihedraltypes=2
new_impropertypes=1
new_types=[2,2,5,2,1]
newmasses=['12.011','15.999']
newpaircoeffs=['0.07000000000000 3.56359487256136 0.07000000000000 3.56359487256136','0.12000000000000 3.02905564167715 0.12000000000000 3.02905564167715']
newbondcoeffs=['525.000 1.2600','200.000 1.4800']
newanglecoeffs=['40.00000 114.00000 50.00000 2.38800','100.00000 132.00000 70.00000 2.22500','50.00000 109.50000 0.00000 0.00000','45.00000 103.00000 0.00000 0.00000','52.00000 108.00000 0.00000 0.00000']
newdihedralcoeffs=['0.0500 6 180 1.00','0.6400 2 180 1.00']
newimpropercoeffs=['96.000 0.00']
new_coeff=[newmasses,newpaircoeffs,newbondcoeffs,newanglecoeffs,newdihedralcoeffs,newimpropercoeffs]
coeffindex=[];
index=['Pair Coeffs','Bond Coeffs','Angle Coeffs','Dihedral Coeffs','Improper Coeffs','Atoms']
#change types
types=data[3:13:2]
for i in range(5):
s=types[i].split()
s[0]=str(int(s[0])+new_types[i])
types[i]=' '.join(s)+'\n'
data[3:13:2]=types
#add coeffs
for i in range(6):
coeffindex.append(findcoeff(data,index[i]))
if re.match('\d.',data[coeffindex[i]-2]):
data=addcoeffs(data,coeffindex[i]-1,new_coeff[i],getcoeffindex(data[coeffindex[i]-2]))
else:
data.insert(coeffindex[i]-1,coeffformat(1,new_coeff[i][0]))
data1 = open("data_new.data", "w")
data1.write("".join(data))
data1.close()
f.close()
The second step can be implemented using LAMMPS’s modeling capabilities. First, note that the models of cellulose and oxidized cellulose are shown below: 
It can be seen that the oxidized form differs from the original only by the absence of the hydroxyl hydrogen atom on the hydroxymethyl group. Interestingly, the two oxygen atoms on the carboxyl carbon after oxidation share the same type, which coincidentally mirrors the two identical hydrogen atoms on the unoxidized hydroxymethyl carbon. Therefore, cellulose oxidation can be achieved simply by deleting the hydroxyl hydrogen and updating the topology information after oxidation (without adding other atoms or redundant topology information, except for adding one out-of-plane dihedral). The topology information before and after oxidation is summarized as follows: Comparison of force field parameters before and after oxidation The different force field parameters were obtained by comparing the polysaccharide data files before and after oxidation generated by CHARMM-GUI. Using LAMMPS’s built-in commands to modify the topological structure of the original data file and assign force field parameters, the CHARMM force field model of the sugar unit after oxidation of hydroxymethyl to carboxyl can be achieved.
units real
atom_style full
dimension 3
boundary f f f
#Force fields
bond_style harmonic
angle_style charmm
dihedral_style charmmfsw
improper_style harmonic
pair_style lj/charmmfsw/coul/charmmfsh 10.0 12.0
pair_modify mix arithmetic
read_data withoutcoeff.data extra/improper/per/atom 5
special_bonds charmm
neighbor 2.0 bin
neigh_modify every 100 delay 0 check yes
include Oxide.coeff
# define groups
group HO id 48 47
group OO id 45 46
group C id 44
#add coeff
# the coeff must be added to data file not the coeff file
# it has been done with python code
#1 delete_atoms
delete_atoms group HO bond yes
#2 set charge
set group OO charge -0.760
set group C charge 0.520
# # # #3 set pair
set group OO type 16
set group C type 15
#4 set bond angle dihedral
########### Bond ###########
group bondOC id 44 45 46
set group bondOC bond 20
group bondCC id 44 29
set group bondCC bond 21
########### Angle ###########
group angle1 id 29 44 45
set group angle1 angle 40
group angle2 id 29 44 46
set group angle2 angle 40
group angle3 id 45 44 46
set group angle3 angle 41
group angle4 id 30 29 44
set group angle4 angle 42
group angle5 id 31 29 44
set group angle5 angle 43
group angle6 id 40 29 44
set group angle6 angle 44
########### Dihedral ###########
group dihedral1 id 30 29 44 45
set group dihedral1 dihedral 91
group dihedral2 id 30 29 44 46
set group dihedral2 dihedral 91
group dihedral3 id 31 29 44 45
set group dihedral3 dihedral 92
group dihedral4 id 31 29 44 46
set group dihedral4 dihedral 92
group dihedral5 id 40 29 44 45
set group dihedral5 dihedral 91
group dihedral6 id 40 29 44 46
set group dihedral6 dihedral 91
#create and set improper
create_bonds single/improper 1 44 29 46 45
write_coeff Oxide.coeff
write_data withoutcoeff.data nocoeff
At this point, the oxidation from hydroxymethyl to carboxyl is complete. The rest is essentially done. To achieve oxidation of the entire cellulose crystal, simply identify the atom indices of the hydroxymethyl carbons in OVITO, write the indices of the atoms that undergo changes as expressions of the selected carbon indices, and randomly select them for oxidation. (A single loop in LAMMPS can handle this.)
Here is an illustration for a feel of it, haha. The principle is the same. Important notes during oxidation:
1 The random number in a variable changes with each call, so create a new variable to store it!!
2 Delete groups promptly!!
units real
atom_style full
dimension 3
boundary f f f
#Force fields
bond_style harmonic
angle_style charmm
dihedral_style charmmfsw
improper_style harmonic
pair_style lj/charmmfsw/coul/charmmfsh 10.0 12.0
pair_modify mix arithmetic
read_data crystal_nocoeff.data extra/improper/per/atom 500
special_bonds charmm
neighbor 2.0 bin
neigh_modify every 100 delay 0 check yes
# write_data crystal_new.data
include crystal.coeff
variable NS loop 30
label loop
variable seed equal ${NS}^2+2*${NS}+3932
variable a290 equal round(random(0,1,${seed})*180)*42+34443
variable a29 equal ${a290}
print " ${a29} "
# variable a29 equal 36921
variable a44 equal ${a29}+14
print " ${a29} "
variable a45 equal ${a29}+15
print " ${a29} "
variable a46 equal ${a29}+16
variable a47 equal ${a29}+17
variable a48 equal ${a29}+18
# define groups
group HO id ${a47} ${a48}
group OO id ${a45} ${a46}
group C id ${a44}
#add coeff
# the coeff must be added to data file not the coeff file
# it has been done with python code
#1 delete_atoms
delete_atoms group HO bond yes
#2 set charge
set group OO charge -0.760
set group C charge 0.520
# # # #3 set pair
set group OO type 12
set group C type 11
#4 set bond angle dihedral
########### Bond ###########
group bondOC id ${a44} ${a45} ${a46}
set group bondOC bond 16
group bondCC id ${a44} ${a29}
set group bondCC bond 17
########### Angle ###########
group angle1 id ${a29} ${a44} ${a45}
set group angle1 angle 33
group angle2 id ${a29} ${a44} ${a46}
set group angle2 angle 33
group angle3 id ${a45} ${a44} ${a46}
set group angle3 angle 34
group angle4 id 30 ${a29} ${a44}
set group angle4 angle 35
group angle5 id 31 ${a29} ${a44}
set group angle5 angle 36
group angle6 id 40 ${a29} ${a44}
set group angle6 angle 37
########### Dihedral ###########
group dihedral1 id 30 ${a29} ${a44} ${a45}
set group dihedral1 dihedral 84
group dihedral2 id 30 ${a29} ${a44} ${a46}
set group dihedral2 dihedral 84
group dihedral3 id 31 ${a29} ${a44} ${a45}
set group dihedral3 dihedral 85
group dihedral4 id 31 ${a29} ${a44} ${a46}
set group dihedral4 dihedral 85
group dihedral5 id 40 ${a29} ${a44} ${a45}
set group dihedral5 dihedral 84
group dihedral6 id 40 ${a29} ${a44} ${a46}
set group dihedral6 dihedral 84
#create and set improper
create_bonds single/improper 1 ${a44} ${a29} ${a46} ${a45}
group HO delete
group OO delete
group C delete
group bondOC delete
group bondCC delete
group angle1 delete
group angle2 delete
group angle3 delete
group angle4 delete
group angle5 delete
group angle6 delete
group dihedral1 delete
group dihedral2 delete
group dihedral3 delete
group dihedral4 delete
group dihedral5 delete
group dihedral6 delete
# delete_atoms group HO bond yes
next NS
jump SELF loop
write_data Oxide.data nocoeff
