Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F109988308
micelle2d.f
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Apr 24, 06:22
Size
5 KB
Mime Type
text/plain
Expires
Sat, Apr 26, 06:22 (1 d, 17 h)
Engine
blob
Format
Raw Data
Handle
25553353
Attached To
rLAMMPS lammps
micelle2d.f
View Options
c Create LAMMPS 2003 input file for 2d LJ simulation of micelles
c
c Syntax: micelle2d < def.micelle2d > data.file
c
c def file contains size of system and number to turn into surfactants
c attaches a random surfactant tail(s) to random head
c if nonflag is set, will attach 2nd-neighbor bonds in surfactant
c solvent atoms = type 1
c micelle heads = type 2
c micelle tails = type 3,4,5,etc.
program micelle2d
parameter (maxatom=100000,maxbond=10000)
real*4 x(2,maxatom)
integer type(maxatom),molecule(maxatom)
integer bondatom(2,maxbond),bondtype(maxbond)
common xprd,yprd,zprd,xboundlo,xboundhi,yboundlo,yboundhi,
$ zboundlo,zboundhi
999 format (3i7,3f8.3)
998 format (4i7)
read (5,*)
read (5,*)
read (5,*) rhostar
read (5,*) iseed
read (5,*) nx,ny
read (5,*) nsurf
read (5,*) r0
read (5,*) ntails
read (5,*) nonflag
natoms = nx*ny
if (natoms+nsurf*ntails.gt.maxatom) then
write (6,*) 'Too many atoms - boost maxatom'
call exit(0)
endif
nbonds = nsurf*ntails
if (nonflag.eq.1) nbonds = nbonds + nsurf*(ntails-1)
if (nbonds.gt.maxbond) then
write (6,*) 'Too many surfactants - boost maxbond'
call exit(0)
endif
c box size
rlattice = (1.0/rhostar) ** 0.5
xboundlo = 0.0
xboundhi = nx*rlattice
yboundlo = 0.0
yboundhi = ny*rlattice
zboundlo = -0.1
zboundhi = 0.1
sigma = 1.0
xprd = xboundhi - xboundlo
yprd = yboundhi - yboundlo
zprd = zboundhi - zboundlo
c initial square lattice of solvents
m = 0
do j = 1,ny
do i = 1,nx
m = m + 1
x(1,m) = xboundlo + (i-1)*rlattice
x(2,m) = yboundlo + (j-1)*rlattice
molecule(m) = 0
type(m) = 1
enddo
enddo
c turn some into surfactants with molecule ID
c head changes to type 2
c create ntails for each head of types 3,4,...
c each tail is at distance r0 away in straight line with random orientation
do i = 1,nsurf
10 m = random(iseed)*natoms + 1
if (m.gt.natoms) m = natoms
if (molecule(m) .ne. 0) goto 10
molecule(m) = i
type(m) = 2
angle = random(iseed)*2.0*3.1415926
do j = 1,ntails
k = (i-1)*ntails + j
x(1,natoms+k) = x(1,m) + cos(angle)*j*r0*sigma
x(2,natoms+k) = x(2,m) + sin(angle)*j*r0*sigma
molecule(natoms+k) = i
type(natoms+k) = 2+j
call pbc(x(1,natoms+k),x(2,natoms+k))
if (j.eq.1) bondatom(1,k) = m
if (j.ne.1) bondatom(1,k) = natoms+k-1
bondatom(2,k) = natoms+k
bondtype(k) = 1
enddo
enddo
c if nonflag is set, add (ntails-1) 2nd nearest neighbor bonds to end
c of bond list
c k = location in bondatom list where nearest neighbor bonds for
c this surfactant are stored
if (nonflag.eq.1) then
nbonds = nsurf*ntails
do i = 1,nsurf
do j = 1,ntails-1
k = (i-1)*ntails + j
nbonds = nbonds + 1
bondatom(1,nbonds) = bondatom(1,k)
bondatom(2,nbonds) = bondatom(2,k+1)
bondtype(nbonds) = 2
enddo
enddo
endif
c write LAMMPS data file
natoms = natoms + nsurf*ntails
nbonds = nsurf*ntails
if (nonflag.eq.1) nbonds = nbonds + nsurf*(ntails-1)
ntypes = 2 + ntails
nbondtypes = 1
if (nonflag.eq.1) nbondtypes = 2
if (nsurf.eq.0) then
ntypes = 1
nbondtypes = 0
endif
write (6,*) 'LAMMPS 2d micelle data file'
write (6,*)
write (6,*) natoms,' atoms'
write (6,*) nbonds,' bonds'
write (6,*) 0,' angles'
write (6,*) 0,' dihedrals'
write (6,*) 0,' impropers'
write (6,*)
write (6,*) ntypes,' atom types'
write (6,*) nbondtypes,' bond types'
write (6,*) 0,' angle types'
write (6,*) 0,' dihedral types'
write (6,*) 0,' improper types'
write (6,*)
write (6,*) xboundlo,xboundhi,' xlo xhi'
write (6,*) yboundlo,yboundhi,' ylo yhi'
write (6,*) zboundlo,zboundhi,' zlo zhi'
write (6,*)
write (6,*) 'Masses'
write (6,*)
do i = 1,ntypes
write (6,*) i,1.0
enddo
write (6,*)
write (6,*) 'Atoms'
write (6,*)
do i = 1,natoms
write (6,999) i,molecule(i),type(i),x(1,i),x(2,i),0.0
enddo
if (nsurf.gt.0) then
write (6,*)
write (6,*) 'Bonds'
write (6,*)
do i = 1,nbonds
write (6,998) i,bondtype(i),bondatom(1,i),bondatom(2,i)
enddo
endif
c write Xmovie bond geometry file
open(1,file='bond.micelle')
write (1,*) 'ITEM: BONDS'
do i = 1,nbonds
write (1,*) bondtype(i),bondatom(1,i),bondatom(2,i)
enddo
close(1)
end
c ************
c Subroutines
c ************
c periodic boundary conditions - map atom back into periodic box
subroutine pbc(x,y)
common xprd,yprd,zprd,xboundlo,xboundhi,yboundlo,yboundhi,
$ zboundlo,zboundhi
if (x.lt.xboundlo) x = x + xprd
if (x.ge.xboundhi) x = x - xprd
if (y.lt.yboundlo) y = y + yprd
if (y.ge.yboundhi) y = y - yprd
return
end
c RNG - compute in double precision, return single
real*4 function random(iseed)
real*8 aa,mm,sseed
parameter (aa=16807.0D0,mm=2147483647.0D0)
sseed = iseed
sseed = mod(aa*sseed,mm)
random = sseed/mm
iseed = sseed
return
end
Event Timeline
Log In to Comment