convert to Fortran 2003
This commit is contained in:
Axel Kohlmeyer
@ -1,233 +0,0 @@
|
||||
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
|
||||
231
tools/micelle2d.f90
Normal file
231
tools/micelle2d.f90
Normal file
@ -0,0 +1,231 @@
|
||||
! Create LAMMPS data file for 2d LJ simulation of micelles
|
||||
!
|
||||
! Syntax: micelle2d < def.micelle2d > data.file
|
||||
!
|
||||
! def file contains size of system and number to turn into surfactants
|
||||
! attaches a random surfactant tail(s) to random head
|
||||
! if nonflag is set, will attach 2nd-neighbor bonds in surfactant
|
||||
! solvent atoms = type 1
|
||||
! micelle heads = type 2
|
||||
! micelle tails = type 3,4,5,etc.
|
||||
|
||||
MODULE box
|
||||
IMPLICIT NONE
|
||||
PUBLIC
|
||||
REAL(KIND=8) :: xprd,yprd,zprd,xboundlo,xboundhi,yboundlo,yboundhi,zboundlo,zboundhi
|
||||
|
||||
CONTAINS
|
||||
|
||||
! periodic boundary conditions - map atom back into periodic box
|
||||
|
||||
SUBROUTINE pbc(x,y)
|
||||
REAL(KIND=8), INTENT(inout) :: x,y
|
||||
|
||||
IF (x < xboundlo) x = x + xprd
|
||||
IF (x >= xboundhi) x = x - xprd
|
||||
IF (y < yboundlo) y = y + yprd
|
||||
IF (y >= yboundhi) y = y - yprd
|
||||
|
||||
END SUBROUTINE pbc
|
||||
END MODULE box
|
||||
|
||||
MODULE rng
|
||||
IMPLICIT NONE
|
||||
|
||||
CONTAINS
|
||||
|
||||
! *very* minimal random number generator
|
||||
|
||||
REAL(KIND=8) FUNCTION random(iseed)
|
||||
IMPLICIT NONE
|
||||
INTEGER, INTENT(inout) :: iseed
|
||||
REAL(KIND=8), PARAMETER :: aa=16807.0_8, mm=2147483647.0_8
|
||||
REAL(KIND=8) :: sseed
|
||||
|
||||
sseed = REAL(iseed)
|
||||
sseed = MOD(aa*sseed,mm)
|
||||
random = sseed/mm
|
||||
iseed = INT(sseed)
|
||||
END FUNCTION random
|
||||
END MODULE rng
|
||||
|
||||
PROGRAM micelle2d
|
||||
USE box
|
||||
USE rng
|
||||
IMPLICIT NONE
|
||||
|
||||
REAL(kind=8), ALLOCATABLE :: x(:,:)
|
||||
INTEGER, ALLOCATABLE :: atomtype(:), molecule(:)
|
||||
INTEGER, ALLOCATABLE :: bondatom(:,:),bondtype(:)
|
||||
INTEGER :: natoms, maxatom, ntypes, nbonds, nbondtypes, iseed
|
||||
INTEGER :: i, j, k, m, nx, ny, nsurf, ntails, nonflag
|
||||
REAL(kind=8) :: rhostar, rlattice, sigma, angle,r0
|
||||
REAL(kind=8), parameter :: pi = 3.14159265358979323846_8
|
||||
LOGICAL :: again
|
||||
|
||||
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
|
||||
maxatom = natoms + nsurf*ntails
|
||||
ALLOCATE(x(2,maxatom), molecule(maxatom), atomtype(maxatom))
|
||||
|
||||
nbonds = nsurf*ntails
|
||||
IF (nonflag.EQ.1) nbonds = nbonds + nsurf*(ntails-1)
|
||||
ALLOCATE(bondatom(2,nbonds), bondtype(nbonds))
|
||||
|
||||
! box size
|
||||
|
||||
rlattice = (1.0_8/rhostar) ** 0.5_8
|
||||
|
||||
xboundlo = 0.0_8
|
||||
xboundhi = nx*rlattice
|
||||
yboundlo = 0.0_8
|
||||
yboundhi = ny*rlattice
|
||||
zboundlo = -0.1_8
|
||||
zboundhi = 0.1_8
|
||||
|
||||
sigma = 1.0_8
|
||||
|
||||
xprd = xboundhi - xboundlo
|
||||
yprd = yboundhi - yboundlo
|
||||
zprd = zboundhi - zboundlo
|
||||
|
||||
! 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
|
||||
atomtype(m) = 1
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
! turn some into surfactants with molecule ID
|
||||
! head changes to type 2
|
||||
! create ntails for each head of types 3,4,...
|
||||
! each tail is at distance r0 away in straight line with random orientation
|
||||
|
||||
DO i = 1,nsurf
|
||||
|
||||
again = .TRUE.
|
||||
DO WHILE(again)
|
||||
m = INT(random(iseed)*natoms + 1)
|
||||
IF (m > natoms) m = natoms
|
||||
IF (molecule(m) /= 0) CYCLE
|
||||
again = .FALSE.
|
||||
END DO
|
||||
molecule(m) = i
|
||||
atomtype(m) = 2
|
||||
|
||||
angle = random(iseed)*2.0_8*pi
|
||||
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
|
||||
atomtype(natoms+k) = 2+j
|
||||
CALL pbc(x(1,natoms+k),x(2,natoms+k))
|
||||
IF (j == 1) bondatom(1,k) = m
|
||||
IF (j /= 1) bondatom(1,k) = natoms+k-1
|
||||
bondatom(2,k) = natoms+k
|
||||
bondtype(k) = 1
|
||||
ENDDO
|
||||
|
||||
ENDDO
|
||||
|
||||
! if nonflag is set, add (ntails-1) 2nd nearest neighbor bonds to end
|
||||
! of bond list
|
||||
! k = location in bondatom list where nearest neighbor bonds for
|
||||
! this surfactant are stored
|
||||
|
||||
IF (nonflag == 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
|
||||
|
||||
! write LAMMPS data file
|
||||
|
||||
natoms = natoms + nsurf*ntails
|
||||
nbonds = nsurf*ntails
|
||||
IF (nonflag == 1) nbonds = nbonds + nsurf*(ntails-1)
|
||||
ntypes = 2 + ntails
|
||||
nbondtypes = 1
|
||||
IF (nonflag == 1) nbondtypes = 2
|
||||
|
||||
IF (nsurf == 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 # molecular'
|
||||
WRITE (6,*)
|
||||
|
||||
DO i = 1,natoms
|
||||
WRITE (6,'(3I7,3F8.3)') i,molecule(i),atomtype(i),x(1,i),x(2,i),0.0
|
||||
ENDDO
|
||||
|
||||
IF (nsurf > 0) THEN
|
||||
|
||||
WRITE (6,*)
|
||||
WRITE (6,*) 'Bonds'
|
||||
WRITE (6,*)
|
||||
|
||||
DO i = 1,nbonds
|
||||
WRITE (6,'(4I7)') i,bondtype(i),bondatom(1,i),bondatom(2,i)
|
||||
ENDDO
|
||||
|
||||
ENDIF
|
||||
|
||||
DEALLOCATE(x,molecule,atomtype,bondtype,bondatom)
|
||||
END PROGRAM micelle2d
|
||||
Reference in New Issue
Block a user