1! Create LAMMPS data file for 2d LJ simulation of micelles
2!
3! Syntax: micelle2d < def.micelle2d > data.file
4!
5! def file contains size of system and number to turn into surfactants
6! attaches a random surfactant tail(s) to random head
7! if nonflag is set, will attach 2nd-neighbor bonds in surfactant
8! solvent atoms = type 1
9! micelle heads = type 2
10! micelle tails = type 3,4,5,etc.
11
12MODULE boxmicelle
13  IMPLICIT NONE
14  PUBLIC
15  REAL(KIND=8) :: xprd,yprd,zprd,xboundlo,xboundhi,yboundlo,yboundhi,zboundlo,zboundhi
16
17CONTAINS
18
19  ! periodic boundary conditions - map atom back into periodic box
20
21  SUBROUTINE pbc(x,y)
22    REAL(KIND=8), INTENT(inout) :: x,y
23
24    IF (x < xboundlo) x = x + xprd
25    IF (x >= xboundhi) x = x - xprd
26    IF (y < yboundlo) y = y + yprd
27    IF (y >= yboundhi) y = y - yprd
28
29  END SUBROUTINE pbc
30END MODULE boxmicelle
31
32MODULE rngmicelle
33  IMPLICIT NONE
34
35CONTAINS
36
37  ! *very* minimal random number generator
38
39  REAL(KIND=8) FUNCTION random(iseed)
40    IMPLICIT NONE
41    INTEGER, INTENT(inout) :: iseed
42    REAL(KIND=8), PARAMETER :: aa=16807.0_8, mm=2147483647.0_8
43    REAL(KIND=8) :: sseed
44
45    sseed = REAL(iseed)
46    sseed = MOD(aa*sseed,mm)
47    random = sseed/mm
48    iseed = INT(sseed)
49  END FUNCTION random
50END MODULE rngmicelle
51
52PROGRAM micelle2d
53  USE boxmicelle
54  USE rngmicelle
55  IMPLICIT NONE
56
57  REAL(kind=8), ALLOCATABLE :: x(:,:)
58  INTEGER, ALLOCATABLE :: atomtype(:), molecule(:)
59  INTEGER, ALLOCATABLE :: bondatom(:,:),bondtype(:)
60  INTEGER :: natoms, maxatom, ntypes, nbonds, nbondtypes, iseed
61  INTEGER :: i, j, k, m, nx, ny, nsurf, ntails, nonflag
62  REAL(kind=8) :: rhostar, rlattice, sigma, angle,r0
63  REAL(kind=8), parameter :: pi = 3.14159265358979323846_8
64  LOGICAL :: again
65
66  READ(5,*)
67  READ(5,*)
68  READ(5,*) rhostar
69  READ(5,*) iseed
70  READ(5,*) nx,ny
71  READ(5,*) nsurf
72  READ(5,*) r0
73  READ(5,*) ntails
74  READ(5,*) nonflag
75
76  natoms = nx*ny
77  maxatom = natoms + nsurf*ntails
78  ALLOCATE(x(2,maxatom), molecule(maxatom), atomtype(maxatom))
79
80  nbonds = nsurf*ntails
81  IF (nonflag.EQ.1) nbonds = nbonds + nsurf*(ntails-1)
82  ALLOCATE(bondatom(2,nbonds), bondtype(nbonds))
83
84! box size
85
86  rlattice = (1.0_8/rhostar) ** 0.5_8
87
88  xboundlo = 0.0_8
89  xboundhi = nx*rlattice
90  yboundlo = 0.0_8
91  yboundhi = ny*rlattice
92  zboundlo = -0.1_8
93  zboundhi = 0.1_8
94
95  sigma = 1.0_8
96
97  xprd = xboundhi - xboundlo
98  yprd = yboundhi - yboundlo
99  zprd = zboundhi - zboundlo
100
101! initial square lattice of solvents
102
103  m = 0
104  DO j = 1,ny
105      DO i = 1,nx
106          m = m + 1
107          x(1,m) = xboundlo + (i-1)*rlattice
108          x(2,m) = yboundlo + (j-1)*rlattice
109          molecule(m) = 0
110          atomtype(m) = 1
111      ENDDO
112  ENDDO
113
114! turn some into surfactants with molecule ID
115!  head changes to type 2
116!  create ntails for each head of types 3,4,...
117!  each tail is at distance r0 away in straight line with random orientation
118
119  DO i = 1,nsurf
120
121      again = .TRUE.
122      DO WHILE(again)
123          m = INT(random(iseed)*natoms + 1)
124          IF (m > natoms) m = natoms
125          IF (molecule(m) /= 0) CYCLE
126          again = .FALSE.
127      END DO
128      molecule(m) = i
129      atomtype(m) = 2
130
131      angle = random(iseed)*2.0_8*pi
132      DO j = 1,ntails
133          k = (i-1)*ntails + j
134          x(1,natoms+k) = x(1,m) + COS(angle)*j*r0*sigma
135          x(2,natoms+k) = x(2,m) + SIN(angle)*j*r0*sigma
136          molecule(natoms+k) = i
137          atomtype(natoms+k) = 2+j
138          CALL pbc(x(1,natoms+k),x(2,natoms+k))
139          IF (j == 1) bondatom(1,k) = m
140          IF (j /= 1) bondatom(1,k) = natoms+k-1
141          bondatom(2,k) = natoms+k
142          bondtype(k) = 1
143      ENDDO
144
145  ENDDO
146
147! if nonflag is set, add (ntails-1) 2nd nearest neighbor bonds to end
148!   of bond list
149! k = location in bondatom list where nearest neighbor bonds for
150!     this surfactant are stored
151
152  IF (nonflag == 1) THEN
153
154      nbonds = nsurf*ntails
155      DO i = 1,nsurf
156          DO j = 1,ntails-1
157              k = (i-1)*ntails + j
158              nbonds = nbonds + 1
159              bondatom(1,nbonds) = bondatom(1,k)
160              bondatom(2,nbonds) = bondatom(2,k+1)
161              bondtype(nbonds) = 2
162          ENDDO
163      ENDDO
164
165  ENDIF
166
167! write LAMMPS data file
168
169  natoms = natoms + nsurf*ntails
170  nbonds = nsurf*ntails
171  IF (nonflag == 1) nbonds = nbonds + nsurf*(ntails-1)
172  ntypes = 2 + ntails
173  nbondtypes = 1
174  IF (nonflag == 1) nbondtypes = 2
175
176  IF (nsurf == 0) THEN
177      ntypes = 1
178      nbondtypes = 0
179  ENDIF
180
181  WRITE (6,*) 'LAMMPS 2d micelle data file'
182  WRITE (6,*)
183
184  WRITE (6,*) natoms,' atoms'
185  WRITE (6,*) nbonds,' bonds'
186  WRITE (6,*) 0,' angles'
187  WRITE (6,*) 0,' dihedrals'
188  WRITE (6,*) 0,' impropers'
189  WRITE (6,*)
190
191  WRITE (6,*) ntypes,' atom types'
192  WRITE (6,*) nbondtypes,' bond types'
193  WRITE (6,*) 0,' angle types'
194  WRITE (6,*) 0,' dihedral types'
195  WRITE (6,*) 0,' improper types'
196  WRITE (6,*)
197
198  WRITE (6,*) xboundlo,xboundhi,' xlo xhi'
199  WRITE (6,*) yboundlo,yboundhi,' ylo yhi'
200  WRITE (6,*) zboundlo,zboundhi,' zlo zhi'
201
202  WRITE (6,*)
203  WRITE (6,*) 'Masses'
204  WRITE (6,*)
205
206  DO i = 1,ntypes
207      WRITE (6,*) i,1.0
208  ENDDO
209
210  WRITE (6,*)
211  WRITE (6,*) 'Atoms # molecular'
212  WRITE (6,*)
213
214  DO i = 1,natoms
215      WRITE (6,'(3I7,3F8.3)') i,molecule(i),atomtype(i),x(1,i),x(2,i),0.0
216  ENDDO
217
218  IF (nsurf > 0) THEN
219
220      WRITE (6,*)
221      WRITE (6,*) 'Bonds'
222      WRITE (6,*)
223
224      DO i = 1,nbonds
225          WRITE (6,'(4I7)') i,bondtype(i),bondatom(1,i),bondatom(2,i)
226      ENDDO
227
228  ENDIF
229
230  DEALLOCATE(x,molecule,atomtype,bondtype,bondatom)
231END PROGRAM micelle2d
232