[gmx-users] building up a long polymer chain
Moeed
lecielll at googlemail.com
Fri Jun 18 21:33:12 CEST 2010
Dear experts,
For hexane molecule I used atom types
opls_157 12.01100 ; all-atom C: CH3 & CH2, alcohols
opls_158 12.01100 ; all-atom C: CH, alcohols
in residue file:
[ HEX ]
[ atoms ]
C1 opls_157 -0.180 1
H11 opls_140 0.060 1
H12 opls_140 0.060 1
H13 opls_140 0.060 1
C2 opls_158 -0.120 2
H21 opls_140 0.060 2
H22 opls_140 0.060 2
C3 opls_158 -0.120 3
H31 opls_140 0.060 3
H32 opls_140 0.060 3
C4 opls_158 -0.120 4
H41 opls_140 0.060 4
H42 opls_140 0.060 4
C5 opls_158 -0.120 5
H51 opls_140 0.060 5
H52 opls_140 0.060 5
C6 opls_157 -0.180 6
H61 opls_140 0.060 6
H62 opls_140 0.060 6
H63 opls_140 0.060 6
now that I am trying to create top file for polyethylene (PE) I doubt if I
pickec right atom types for hexane since I see in ffoplsaa.atp file:
opls_135 12.01100 ; alkane CH3
opls_136 12.01100 ; alkane CH2
for alkanes. Please let me know if I am right and I should use
opls_157 12.01100 ; all-atom C: CH3 & CH2, alcohols
opls_158 12.01100 ; all-atom C: CH, alcohols
for PE. or maybe this dies not affect the results since charges and mass are
the same..?!
Also: can opls_157 be used for both CH3 & CH2,on PE chain? what is the
suitable atom type for carbon atoms on PE chain?
2- I spent alot of time finding a solution for the following problem: I am
trying to build up a very long PE chain. to do this I am using the following
residues in rtp file:
; Polyethylene - this is an internal residue
[ Eth ]
[ atoms ]
C1 opls_136 -0.120 1
H11 opls_140 0.060 1
H12 opls_140 0.060 1
C2 opls_136 -0.120 2
H21 opls_140 0.060 2
H22 opls_140 0.060 2
[ bonds ]
C1 -C2
C1 H11
C1 H12
C1 C2
C2 H21
C2 H22
C2 +C1
; Terminal PE residue ("beginning" of chain)
; designation arbitrary, C1 is -CH3
; designation arbitrary, C1 is -CH3
[ EthB ]
[ atoms ]
C1 opls_135 -0.180 1
H11 opls_140 0.060 1
H12 opls_140 0.060 1
H13 opls_140 0.060 1
C2 opls_136 -0.120 2
H21 opls_140 0.060 2
H22 opls_140 0.060 2
[ bonds ]
C1 H11
C1 H12
C1 H13
C1 C2
C2 H21
C2 H22
C2 +C
; Terminal PE residue ("end" of chain)
; designation arbitrary, C2 is -CH3
[ EthE ]
[ atoms ]
C1 opls_136 -0.120 1
H11 opls_140 0.060 1
H12 opls_140 0.060 1
C2 opls_135 -0.180 2
H21 opls_140 0.060 2
H22 opls_140 0.060 2
H23 opls_140 0.060 2
[ bonds ]
C1 -C
C1 H11
C1 H12
C1 C2
C2 H21
C2 H22
C2 H23
However, the software I am using(Ascalaph Desingner) to build up PE chain
lists atoms as shown below: (I changed the res names to Eth and EthB
myself). i.e. atoms have no numbering meaning that if I want to use the
above resudies in rtp file I will have to go thorough the whole structure
file(below) and do the numbering by myself (H to H11, H12, H21, H22...)and C
to C1 C2) which is a tedious work. when I want to generate a molecule of
40000 Mw (about 1400 monomer units! or about 8000 atoms)
HETATM 1 CB EthB 1 2.739 2.554 -0.012 0.00 -0.39
C
HETATM 2 H EthB 1 2.736 3.644 -0.012 0.00 0.13
H
HETATM 3 H EthB 1 2.190 2.217 -0.893 0.00 0.13
H
HETATM 4 H EthB 1 2.190 2.217 0.868 0.00 0.13
H
HETATM 5 C Eth 2 4.155 1.993 -0.012 0.00 -0.26
C
HETATM 6 H Eth 2 4.696 2.362 -0.887 0.00 0.13
H
HETATM 7 H Eth 2 4.696 2.362 0.862 0.00 0.13
H
HETATM 8 C Eth 3 4.163 0.470 -0.012 0.00 -0.26
C
HETATM 9 H Eth 3 3.623 0.101 0.862 0.00 0.13
H
HETATM 10 H Eth 3 3.623 0.101 -0.887 0.00 0.13
H
HETATM 11 C Eth 4 5.579 -0.091 -0.012 0.00 -0.26
C
HETATM 12 H Eth 4 6.119 0.278 -0.887 0.00 0.13
H
.
.
I numbered 350 atoms (60 monomer units) and simulation is working now.
56
1EthB C1 1 0.860 2.362 1.500
1EthB H11 2 0.860 2.471 1.500
1EthB H12 3 0.805 2.329 1.412
1EthB H13 4 0.805 2.329 1.588
1EthB C2 5 1.002 2.306 1.500
1EthB H21 6 1.056 2.343 1.412
1EthB H22 7 1.056 2.343 1.587
2Eth C1 8 1.002 2.154 1.500
2Eth H11 9 0.948 2.117 1.587
2Eth H12 10 0.948 2.117 1.412
2Eth C2 11 1.144 2.098 1.500
2Eth H21 12 1.198 2.135 1.412
2Eth H22 13 1.198 2.135 1.587
3Eth C1 14 1.145 1.946 1.500
3Eth H11 15 1.091 1.909 1.587
3Eth H12 16 1.091 1.909 1.412
3Eth C2 17 1.286 1.890 1.500
3Eth H21 18 1.340 1.926 1.412
3Eth H22 19 1.340 1.926 1.587
4Eth C1 20 1.287 1.737 1.500
4Eth H11 21 1.233 1.700 1.587
4Eth H12 22 1.233 1.700 1.412
4Eth C2 23 1.429 1.681 1.500
4Eth H21 24 1.483 1.718 1.412
4Eth H22 25 1.483 1.718 1.587
5Eth C1 26 1.429 1.529 1.500
Nevertheless, I am trying to find out a better way to perform simulation for
very longer chain. I examined x2top command to het top file with no atoms
numbers (structure file .pdb above) and got gro file using editconf :
56
1EthB CB 1 0.860 2.362 1.500
1EthB H 2 0.859 2.471 1.500
1EthB H 3 0.805 2.328 1.412
1EthB H 4 0.805 2.328 1.588
2Eth C 5 1.001 2.306 1.500
2Eth H 6 1.055 2.343 1.413
2Eth H 7 1.055 2.343 1.587
3Eth C 8 1.002 2.153 1.500
3Eth H 9 0.948 2.116 1.587
3Eth H 10 0.948 2.116 1.413
4Eth C 11 1.144 2.097 1.500
4Eth H 12 1.198 2.134 1.413
4Eth H 13 1.198 2.134 1.587
5Eth C 14 1.145 1.945 1.500
5Eth H 15 1.090 1.908 1.587
5Eth H 16 1.090 1.908 1.413
6Eth C 17 1.286 1.889 1.500
6Eth H 18 1.340 1.926 1.413
6Eth H 19 1.340 1.926 1.587
7Eth C 20 1.287 1.736 1.500
7Eth H 21 1.233 1.700 1.587
The top file I got has contains exactly the same listing of [bonds] [pairs]
[angles] [dihedrals] as the one I got from the structure file for which I
did the numbering by hand. But I am unsure of the forcefield parameters if
they are extracted from itp file properly because after doing energy
minimization step I get LINCS error.
Back Off! I just backed up ener.edr to ./#ener.edr.1#
Steepest Descents:
Tolerance (Fmax) = 1.00000e+03
Number of steps = 200
Step= 0, Dmax= 1.0e-02 nm, Epot= 5.25663e+04 Fmax= 2.91635e+03, atom= 53
Step= 1, Dmax= 1.0e-02 nm, Epot= 5.20487e+04 Fmax= 6.27470e+03, atom= 1
Step= 2, Dmax= 1.2e-02 nm, Epot= 5.12076e+04 Fmax= 5.74320e+03, atom= 56
Step= 3, Dmax= 1.4e-02 nm, Epot= 5.08573e+04 Fmax= 4.77251e+03, atom= 7
Step= 4, Dmax= 1.7e-02 nm, Epot= 5.04166e+04 Fmax= 7.91400e+03, atom= 53
Step= 5, Dmax= 2.1e-02 nm, Epot= 4.88521e+04 Fmax= 1.11208e+04, atom= 7
Step= 6, Dmax= 2.5e-02 nm, Epot= 4.66591e+04 Fmax= 8.34560e+03, atom= 48
Step 7, time 0.014 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.005015, max 0.011148 (between atoms 5 and 6)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
Step= 7, Dmax= 3.0e-02 nm, Epot= 4.41164e+04 Fmax= 1.36155e+04, atom= 23
Step= 8, Dmax= 3.6e-02 nm, Epot= 4.22917e+04 Fmax= 8.05183e+03, atom= 23
Step 9, time 0.018 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.022162, max 0.109777 (between atoms 23 and 25)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
23 25 51.2 0.1086 0.1210 0.1090
23 24 44.4 0.1090 0.1150 0.1090
8 10 32.1 0.1100 0.1113 0.1090
Step= 10, Dmax= 2.1e-02 nm, Epot= 4.13912e+04 Fmax= 1.80418e+04, atom= 23
Step= 11, Dmax= 2.6e-02 nm, Epot= 4.03985e+04 Fmax= 6.68523e+03, atom= 39
Step= 12, Dmax= 3.1e-02 nm, Epot= 3.98609e+04 Fmax= 3.40293e+04, atom= 23
Step= 13, Dmax= 3.7e-02 nm, Epot= 3.88153e+04 Fmax= 6.68366e+03, atom= 23
.
.
.
Q 2-1 Please kindly furnish me with your comments on this issue. *I need a
residue for internal ethylene monomer in rtp file which depends as less as
possible on atoms numbering (H11,...). *I think the solution will help many
gmx users dealing with polymers. what would be the best approach?
my current internal repeating unit has numbering of atoms:
[Eth ]
[ atoms ]
C1 opls_136 -0.120 1
H11 opls_140 0.060 1
H12 opls_140 0.060 1
C2 opls_136 -0.120 2
H21 opls_140 0.060 2
H22 opls_140 0.060 2
[ bonds ]
C1 -C2
C1 H11
C1 H12
C1 C2
C2 H21
C2 H22
C2 +C1
Q 2-2. Does x2top command have advantage over pdb2gmx command with regards
to atom numbering?
Many many thanks,
Best,
Moeed
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20100618/33a50103/attachment.html>
More information about the gromacs.org_gmx-users
mailing list