← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to tools/moltemplate/examples/coarse_grained_examples/vesicle_Brannigan2005+Bellesia2010/moltemplate_files/CGLipidBr2005.lt

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2015-04-29 23:44:49 UTC
  • mfrom: (5.1.3 experimental)
  • Revision ID: package-import@ubuntu.com-20150429234449-mbhy9utku6hp6oq8
Tags: 0~20150313.gitfa668e1-1
Upload into unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
tools/moltemplate/examples/coarse_grained_examples/vesicle_Brannigan2005+Bellesia2010/moltemplate_files/CGLipidBr2005.lt
 
1
# Note:
 
2
#
 
3
#    This example may require additional features to be added to LAMMPS.  If
 
4
# LAMMPS complains about an "Invalid pair_style", then download copy the
 
5
# "additional_lammps_code" from moltemplate.org, unpack it into your LAMMPS
 
6
# "src" directory and recompile LAMMPS.
 
7
#
 
8
# -------- Description --------
 
9
#
 
10
# This example contains an implementation of the DPPC lipid bilayer described in
 
11
#      G. Brannigan, P.F. Philips, and F.L.H. Brown,
 
12
#      Physical Review E, Vol 72, 011915 (2005)
 
13
# and:
 
14
#      M.C. Watson, E.S. Penev, P.M. Welch, and F.L.H. Brown
 
15
#      J. Chem. Phys. 135, 244701 (2011)
 
16
#
 
17
# As in Watson(JCP 2011), rigid bond-length constraints have been replaced 
 
18
# by harmonic bonds.
 
19
#
 
20
#   A truncated version of this lipid (named "DLPC") has also been added.
 
21
# Unlike the original "DPPC" molecule model, "DLPC" has not been carefully
 
22
# parameterized to reproduce the correct behavior in a lipid bilayer mixture.
 
23
#
 
24
#   Units:
 
25
#
 
26
# The "epsilon" parameter in their model is approximately 2.75 kJ/mole
 
27
#     ( = 0.657265774378585 kCal/mole, using 1kCal=4.184kJ)
 
28
# The "sigma" parameter corresponds to 7.5 angstroms.
 
29
 
 
30
 
 
31
CGLipidBr2005 {
 
32
 
 
33
 
 
34
  write_once("In Init") {
 
35
    # -- Default styles for "CGLipidBr2005" --
 
36
    units           real
 
37
    atom_style      full
 
38
    # (Hybrid force field styles were used for portability.)
 
39
    bond_style      hybrid harmonic
 
40
 
 
41
    #angle_style     hybrid cosine/delta # <- used in the original article
 
42
    angle_style     hybrid harmonic  # <- prevents unphysical acute angle turns
 
43
    # Explanation:
 
44
    # angle_style cosine/delta:  U(theta) = k*(1-cos(theta-theta0))
 
45
    # angle_style     harmonic:  U(theta) = k*(theta-theta0)^2
 
46
 
 
47
    dihedral_style  none
 
48
    improper_style  none
 
49
    pair_style  hybrid table linear 1130 &
 
50
                lj/charmm/coul/charmm/inter es4k4l 14.5 15
 
51
    pair_modify     mix arithmetic
 
52
    special_bonds   lj 0.0 0.0 1.0   # turn off pairs if "less than 3 bonds"
 
53
 
 
54
    neighbor 2.0 multi        # <- perhaps unnecessary
 
55
    comm_modify mode multi    # <- perhaps unnecessary
 
56
  }
 
57
 
 
58
 
 
59
  DPPC {
 
60
    write("Data Atoms") {
 
61
      $atom:h  $mol:. @atom:head     0.0  0.00  0.00  33.75 # DPPC head atom
 
62
      $atom:i  $mol:. @atom:../int   0.0 -1.00  0.00  26.25
 
63
      $atom:t1 $mol:. @atom:../tail  0.0  1.00  0.00  18.75
 
64
      $atom:t2 $mol:. @atom:../tail  0.0 -1.00  0.00  11.25
 
65
      $atom:t3 $mol:. @atom:../tail  0.0  1.00  0.00   3.75
 
66
    }
 
67
    write("Data Bonds") {
 
68
      $bond:b1    @bond:../backbone   $atom:h  $atom:i
 
69
      $bond:b2    @bond:../backbone   $atom:i  $atom:t1
 
70
      $bond:b3    @bond:../backbone   $atom:t1 $atom:t2
 
71
      $bond:b4    @bond:../backbone   $atom:t2 $atom:t3
 
72
    }
 
73
    write("Data Angles") {
 
74
      $angle:a1   @angle:../backbone  $atom:h   $atom:i   $atom:t1
 
75
      $angle:a2   @angle:../backbone  $atom:i   $atom:t1  $atom:t2
 
76
      $angle:a3   @angle:../backbone  $atom:t1  $atom:t2  $atom:t3
 
77
    }
 
78
 
 
79
    # Define properties of the local (lipid-specific) atom:head type atom:
 
80
    write_once("Data Masses") {
 
81
      @atom:head  200.0
 
82
    }
 
83
    write_once("In Settings") {
 
84
      pair_coeff @atom:head @atom:head   lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
 
85
      pair_coeff @atom:../int @atom:head lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
 
86
    }
 
87
 
 
88
  } #DPPC
 
89
 
 
90
 
 
91
  DLPC {
 
92
    write("Data Atoms") {
 
93
      $atom:h  $mol:. @atom:head     0.0  0.00  0.00  30.00 # DLPC head atom
 
94
      $atom:i  $mol:. @atom:../int   0.0 -1.00  0.00  22.50
 
95
      $atom:t1 $mol:. @atom:../tail  0.0  1.00  0.00  15.00
 
96
      $atom:t2 $mol:. @atom:../tail  0.0 -1.00  0.00   7.50
 
97
    }
 
98
    write("Data Bonds") {
 
99
      $bond:b1    @bond:../backbone   $atom:h  $atom:i
 
100
      $bond:b2    @bond:../backbone   $atom:i  $atom:t1
 
101
      $bond:b3    @bond:../backbone   $atom:t1 $atom:t2
 
102
    }
 
103
    write("Data Angles") {
 
104
      $angle:a1   @angle:../backbone  $atom:h   $atom:i   $atom:t1
 
105
      $angle:a2   @angle:../backbone  $atom:i   $atom:t1  $atom:t2
 
106
    }
 
107
    # Define properties of the local (lipid-specific) atom:head type atom:
 
108
    write_once("Data Masses") {
 
109
      @atom:head  200.0
 
110
    }
 
111
    write_once("In Settings") {
 
112
      pair_coeff @atom:head @atom:head lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
 
113
      pair_coeff @atom:../int @atom:head lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
 
114
    }
 
115
  } #DLPC
 
116
 
 
117
 
 
118
  # Particles and properties shared by all lipid types:
 
119
 
 
120
  write_once("Data Masses") {
 
121
    @atom:int     200.0
 
122
    @atom:tail    200.0
 
123
    @atom:head    200.0  #<- Default head type.  We may override it later.
 
124
  }
 
125
 
 
126
  write_once("In Settings") {
 
127
    # -- Default settings/parameters for "CGLipidBr2005" --
 
128
    # (Hybrid bond & angle styles were used for portability.)
 
129
 
 
130
    #  As in Watson(JCP 2011), rigid bond-length constraints 
 
131
    #  have been replaced by harmonic bonds.
 
132
    #  The k_theta parameter should lie in between 5*epsilon and 10*epsilon.
 
133
    bond_coeff   @bond:backbone   harmonic     116.847 7.5 #<--2*5000*eps/sig^2
 
134
  }
 
135
 
 
136
  write_once("In Settings") {
 
137
    #  cosine/delta: U(theta) = k*(1-cos(theta-theta0))
 
138
    #angle_coeff  @angle:backbone  cosine/delta  4.60086042   180 #<--  7*eps
 
139
    #  harmonic: U(theta) = k*(theta-theta0)^2     not (k/2)*(theta-theta0)^2
 
140
    angle_coeff  @angle:backbone  harmonic      9.85898661    180 #<-->30*eps
 
141
  }
 
142
  # I use a stiffer bond-angle than the original Brannigan & Brown 2005 paper
 
143
  # to attempt to compensate for the fact that here we are using a lipid
 
144
  # mixture of DPPC and DLPC.  (The mixture of lipids introduces a great deal 
 
145
  # of disorder into the bilayer which would not be present in a DPPC bilayer.
 
146
  # This causes pores to form.  Increasing the angle stiffness prevents this.)
 
147
 
 
148
  write_once("In Settings") {
 
149
 
 
150
    # The interaction of "atom:int" with other "atom:int" atoms is given by
 
151
    # epsilon*(0.4*(sigma/r)^12 - 3.0*(sigma/r)^2), shifted and cutoff at 
 
152
    # r=3*sigma.  This was implemented using pair_style table.
 
153
    # Unfortunately, mixing lj/charmm and "table" pair styles in the same
 
154
    # simulation is very inneficient.
 
155
 
 
156
  pair_coeff @atom:int  @atom:int  table  table_int.dat  INT
 
157
 
 
158
    # The interaction of tail beads with eachother is given by the formula below
 
159
    # and with other atoms ...using Lorenz-Berthelot and "repulsive wins" rules:
 
160
    # epsilon*(0.4*(sigma/r)^12 - 1.0*(sigma/r)^6),
 
161
  pair_coeff @atom:tail @atom:tail lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 -1
 
162
  pair_coeff @atom:int  @atom:tail lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 -1
 
163
 
 
164
    # The interaction of head beads which all other beads is given by:
 
165
    # epsilon*(0.4*(sigma/r)^12 - 0.0*(sigma/r)^6),
 
166
 
 
167
  pair_coeff @atom:head @atom:head lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
 
168
  pair_coeff @atom:int  @atom:head lj/charmm/coul/charmm/inter 0.1643 7.5 0.4 0
 
169
 
 
170
 
 
171
  }  # write_once("In Settings")
 
172
 
 
173
  # Note: I divided epsilon by 4 to get "0.1643" because we are using the 
 
174
  # "es4k4l" coeffstyle, corresponding to U(r)=eps(4*K*(s/r)^12 + 4*L*(s/r)^6)
 
175
  # (The "es4k4l" coeffstyle is the default.) Using this convention makes it
 
176
  # easier to mix this coarse-grained lipid model with other molecular models.
 
177
 
 
178
 
 
179
 
 
180
} # CGLipidBr2005
 
181
 
 
182
 
 
183
 
 
184
 
 
185
 
 
186
 
 
187
 
 
188
 
 
189
# Note: This example has not been optimized for speed.
 
190
#
 
191
# Unfortunately, using both lj/charmm and "table" pair styles in the same
 
192
# simulation seems to be very inneficient.  (The simulation is twice as slow
 
193
# as using only the "lj/charmm" pair styles for every pairwise interaction, 
 
194
# ...and about 25% slower than using "table" for every pairwise interaction.
 
195
# However the lennard-jones pair styles support mixing, so we use them to
 
196
# make it easier to run these molecules with other molecules which don't use
 
197
# pair_table.  I felt that portability was worth the extra 25% slow down.)
 
198