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