1! Create LAMMPS data file for collection of
2!   polymer bead-spring chains of various lengths and bead sizes
3! Syntax: chain < def.chain > data.file
4!   def.chain is input file that specifies the chains
5!   data.file is output file that will be input for LAMMPS
6! includes image flags in data file so chains can be unraveled later
7
8MODULE boxchain
9  IMPLICIT NONE
10  PUBLIC
11  REAL(KIND=8) :: xprd,yprd,zprd,xboundlo,xboundhi,yboundlo,yboundhi,zboundlo,zboundhi
12
13CONTAINS
14
15  ! periodic boundary conditions - map atom back into periodic box
16
17  SUBROUTINE pbc(x,y,z)
18    REAL(KIND=8), INTENT(inout) :: x,y,z
19
20    IF (x < xboundlo) x = x + xprd
21    IF (x >= xboundhi) x = x - xprd
22    IF (y < yboundlo) y = y + yprd
23    IF (y >= yboundhi) y = y - yprd
24    IF (z < zboundlo) z = z + zprd
25    IF (z >= zboundhi) z = z - zprd
26
27  END SUBROUTINE pbc
28END MODULE boxchain
29
30MODULE rngchain
31  IMPLICIT NONE
32
33CONTAINS
34
35  ! *very* minimal random number generator
36
37  REAL(KIND=8) FUNCTION random(iseed)
38    IMPLICIT NONE
39    INTEGER, INTENT(inout) :: iseed
40    REAL(KIND=8), PARAMETER :: aa=16807.0_8, mm=2147483647.0_8
41    REAL(KIND=8) :: sseed
42
43    sseed = REAL(iseed)
44    sseed = MOD(aa*sseed,mm)
45    random = sseed/mm
46    iseed = INT(sseed)
47  END FUNCTION random
48END MODULE rngchain
49
50PROGRAM chain
51  USE boxchain
52  USE rngchain
53  IMPLICIT NONE
54
55  INTEGER, ALLOCATABLE :: nchain(:),nmonomer(:)
56  INTEGER, ALLOCATABLE :: ntype(:),nbondtype(:)
57  INTEGER, ALLOCATABLE :: atomtype(:),molecule(:)
58  INTEGER, ALLOCATABLE :: imagex(:),imagey(:),imagez(:)
59  REAL(KIND=8), ALLOCATABLE :: x(:),y(:),z(:)
60  REAL(KIND=8), ALLOCATABLE :: bondlength(:),restrict(:)
61  INTEGER :: i, n, m, nmolecule, natoms, ntypes, nbonds, nbondtypes
62  INTEGER :: swaptype, iseed, nsets, iset, ichain, imonomer
63  REAL(KIND=8) :: r, rhostar, volume, rsq, xinner, yinner, zinner, xsurf, ysurf, zsurf
64  REAL(KIND=8) :: x0, y0, z0, x1, y1, z1, x2, y2, z2, dx, dy, dz
65
66  LOGICAL :: again
67
68  ! read chain definitions
69
70  READ (5,*)
71  READ (5,*)
72  READ (5,*) rhostar
73  READ (5,*) iseed
74  READ (5,*) nsets
75  READ (5,*) swaptype
76
77  ALLOCATE(nchain(nsets))
78  ALLOCATE(nmonomer(nsets))
79  ALLOCATE(ntype(nsets))
80  ALLOCATE(nbondtype(nsets))
81  ALLOCATE(bondlength(nsets))
82  ALLOCATE(restrict(nsets))
83
84  DO iset = 1,nsets
85      READ (5,*)
86      READ (5,*) nchain(iset)
87      READ (5,*) nmonomer(iset)
88      READ (5,*) ntype(iset)
89      READ (5,*) nbondtype(iset)
90      READ (5,*) bondlength(iset)
91      READ (5,*) restrict(iset)
92  ENDDO
93
94  ! natoms = total # of monomers
95
96  natoms = 0
97  DO iset = 1,nsets
98      natoms = natoms + nchain(iset)*nmonomer(iset)
99  ENDDO
100
101  ALLOCATE(x(natoms))
102  ALLOCATE(y(natoms))
103  ALLOCATE(z(natoms))
104  ALLOCATE(atomtype(natoms))
105  ALLOCATE(molecule(natoms))
106  ALLOCATE(imagex(natoms))
107  ALLOCATE(imagey(natoms))
108  ALLOCATE(imagez(natoms))
109
110  ! setup box size (sigma = 1.0)
111
112  volume = natoms/rhostar
113  xprd = volume**(1.0/3.0)
114  yprd = xprd
115  zprd = xprd
116
117  xboundlo = -xprd/2.0
118  xboundhi = -xboundlo
119  yboundlo = xboundlo
120  yboundhi = xboundhi
121  zboundlo = xboundlo
122  zboundhi = xboundhi
123
124  ! generate random chains
125  ! loop over sets and chains in each set
126
127  n = 0
128  nmolecule = 0
129
130  DO iset = 1,nsets
131      DO ichain = 1,nchain(iset)
132          nmolecule = nmolecule + 1
133
134          ! random starting point for the chain in the box
135
136          x1 = 0.0
137          y1 = 0.0
138          z1 = 0.0
139          x2 = xboundlo + random(iseed)*xprd
140          y2 = yboundlo + random(iseed)*yprd
141          z2 = zboundlo + random(iseed)*zprd
142
143          ! store 1st monomer of chain
144          ! 1st monomer is always in original box (image = 0)
145
146          CALL pbc(x2,y2,z2)
147
148          n = n + 1
149          x(n) = x2
150          y(n) = y2
151          z(n) = z2
152          atomtype(n) = ntype(iset)
153          imagex(n) = 0
154          imagey(n) = 0
155          imagez(n) = 0
156          IF (swaptype == 0) THEN
157              molecule(n) = nmolecule
158          ELSE
159              molecule(n) = 1
160          END IF
161
162          ! generate rest of monomers in this chain
163
164          DO imonomer = 2, nmonomer(iset)
165
166              x0 = x1
167              y0 = y1
168              z0 = z1
169              x1 = x2
170              y1 = y2
171              z1 = z2
172
173              again = .TRUE.
174              DO WHILE (again)
175                  ! random point inside sphere of unit radius
176
177                  xinner = 2.0*random(iseed) - 1.0
178                  yinner = 2.0*random(iseed) - 1.0
179                  zinner = 2.0*random(iseed) - 1.0
180                  rsq = xinner*xinner + yinner*yinner + zinner*zinner
181                  IF (rsq > 1.0) CYCLE
182
183                  ! project point to surface of sphere of unit radius
184
185                  r = SQRT(rsq)
186                  xsurf = xinner/r
187                  ysurf = yinner/r
188                  zsurf = zinner/r
189
190                  ! create new point by scaling unit offsets by bondlength (sigma = 1.0)
191
192                  x2 = x1 + xsurf*bondlength(iset)
193                  y2 = y1 + ysurf*bondlength(iset)
194                  z2 = z1 + zsurf*bondlength(iset)
195
196                  ! check that new point meets restriction requirement
197                  ! only for 3rd monomer and beyond
198
199                  dx = x2 - x0
200                  dy = y2 - y0
201                  dz = z2 - z0
202                  r = SQRT(dx*dx + dy*dy + dz*dz)
203
204                  IF (imonomer > 2 .AND. r <= restrict(iset)) CYCLE
205
206                  ! store new point
207                  again = .FALSE.
208
209                  ! if delta to previous bead is large, then increment/decrement image flag
210
211                  CALL pbc(x2,y2,z2)
212                  n = n + 1
213                  x(n) = x2
214                  y(n) = y2
215                  z(n) = z2
216                  atomtype(n) = ntype(iset)
217
218                  IF (ABS(x(n)-x(n-1)) < 2.0*bondlength(iset)) THEN
219                      imagex(n) = imagex(n-1)
220                  ELSE IF (x(n) - x(n-1) < 0.0) THEN
221                      imagex(n) = imagex(n-1) + 1
222                  ELSE IF (x(n) - x(n-1) > 0.0) THEN
223                      imagex(n) = imagex(n-1) - 1
224                  ENDIF
225
226                  IF (ABS(y(n)-y(n-1)) < 2.0*bondlength(iset)) THEN
227                      imagey(n) = imagey(n-1)
228                  ELSE IF (y(n) - y(n-1) < 0.0) THEN
229                      imagey(n) = imagey(n-1) + 1
230                  ELSE IF (y(n) - y(n-1) > 0.0) THEN
231                      imagey(n) = imagey(n-1) - 1
232                  ENDIF
233
234                  IF (ABS(z(n)-z(n-1)) < 2.0*bondlength(iset)) THEN
235                      imagez(n) = imagez(n-1)
236                  ELSE IF (z(n) - z(n-1) < 0.0) THEN
237                      imagez(n) = imagez(n-1) + 1
238                  ELSE IF (z(n) - z(n-1) > 0.0) THEN
239                      imagez(n) = imagez(n-1) - 1
240                  ENDIF
241
242                  IF (swaptype == 0) THEN
243                      molecule(n) = nmolecule
244                  ELSE IF (swaptype == 1) THEN
245                      molecule(n) = imonomer
246                  ELSE IF (swaptype == 2) THEN
247                      IF (imonomer <= nmonomer(iset)/2) THEN
248                          molecule(n) = imonomer
249                      ELSE
250                          molecule(n) = nmonomer(iset)+1-imonomer
251                      ENDIF
252                  ENDIF
253              ENDDO
254          ENDDO
255
256      ENDDO
257  ENDDO
258
259  ! compute quantities needed for LAMMPS file
260
261  nbonds = 0
262  ntypes = 0
263  nbondtypes = 0
264  DO iset = 1,nsets
265      nbonds = nbonds + nchain(iset)*(nmonomer(iset)-1)
266      IF (ntype(iset) > ntypes) ntypes = ntype(iset)
267      IF (nbondtype(iset) > nbondtypes) nbondtypes = nbondtype(iset)
268  ENDDO
269
270  ! write out LAMMPS file
271
272  WRITE (6,*) 'LAMMPS FENE chain data file'
273  WRITE (6,*)
274
275  WRITE (6,*) natoms,' atoms'
276  WRITE (6,*) nbonds,' bonds'
277  WRITE (6,*) 0,' angles'
278  WRITE (6,*) 0,' dihedrals'
279  WRITE (6,*) 0,' impropers'
280  WRITE (6,*)
281
282  WRITE (6,*) ntypes,' atom types'
283  WRITE (6,*) nbondtypes,' bond types'
284  WRITE (6,*) 0,' angle types'
285  WRITE (6,*) 0,' dihedral types'
286  WRITE (6,*) 0,' improper types'
287  WRITE (6,*)
288
289  WRITE (6,'(2F15.6,A)') xboundlo,xboundhi,' xlo xhi'
290  WRITE (6,'(2F15.6,A)') yboundlo,yboundhi,' ylo yhi'
291  WRITE (6,'(2F15.6,A)') zboundlo,zboundhi,' zlo zhi'
292
293  WRITE (6,*)
294  WRITE (6,*) 'Masses'
295  WRITE (6,*)
296
297  DO i = 1,ntypes
298      WRITE (6,'(i3,f5.1)') i,1.0
299  ENDDO
300
301  WRITE (6,*)
302  WRITE (6,*) 'Atoms # molecular'
303  WRITE (6,*)
304
305  DO i = 1,natoms
306      WRITE (6,'(I10,I8,I8,3F10.4,3I4)') i,molecule(i),atomtype(i),x(i),y(i),z(i), &
307          imagex(i),imagey(i),imagez(i)
308  ENDDO
309
310  IF (nbonds > 0) THEN
311      WRITE (6,*)
312      WRITE (6,*) 'Bonds'
313      WRITE (6,*)
314
315      n = 0
316      m = 0
317      DO iset = 1,nsets
318          DO ichain = 1,nchain(iset)
319              DO imonomer = 1,nmonomer(iset)
320                  n = n + 1
321                  IF (imonomer /= nmonomer(iset)) THEN
322                      m = m + 1
323                      WRITE (6,'(i9,i3,2i9)') m,nbondtype(iset),n,n+1
324                  ENDIF
325              ENDDO
326          ENDDO
327      ENDDO
328  ENDIF
329
330  DEALLOCATE(nchain, nmonomer, ntype, nbondtype, bondlength, restrict)
331  DEALLOCATE(x, y, z, atomtype, molecule, imagex, imagey, imagez)
332
333END PROGRAM chain
334