Homec4science
Diffusion lammps (lammm-master) Browse

lammps/tools/createatoms/createAtoms.f
4df5a5dee043lammm-master

Last Change
Raw File

createAtoms.f
View Options

c X. W. Zhou, xzhou@sandia.gov
include 'createAtoms.h'
call latgen()
call delete()
call hitvalue()
call disturb()
call systemshift()
call colect()
if (ilatseed .ne. 0) call randomize()
call velgen()
call outputfiles()
if (float(nntype(1)) .ne. 0.0)
* write(6,*)(float(nntype(i))/float(nntype(1)),i=1,ntypes)
write(6,*)'atoms of different species ',(nntype(i),i=1,ntypes)
stop
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine latgen()
include 'createAtoms.h'
character*8 lattype
common /latinfo/ lattype,alat(3),xrot(3),yrot(3),zrot(3),
* periodicity(3),xbound(2),ybound(2),zbound(2),strain(3),
* delx(3),rcell(3),ccell(nelmax)
common /iran/ iseed
namelist /maincard/ ntypes,amass,ielement,ilatseed,
* perlb,perub,iseed,xy,xz,yz
namelist /latcard/ lattype,alat,xrot,yrot,zrot,
* periodicity,xbound,ybound,zbound,
* strain,delx
iseed=3245566
natoms=0
ilatseed=0
xy=1.0e12
xz=1.0e12
yz=1.0e12
read(5,maincard)
do i = 1,3
perlen(i)=perub(i)-perlb(i)
enddo
do i=1,ntypes
nntype(i)=0
enddo
1 continue
strain(1)=0.0
strain(2)=0.0
strain(3)=0.0
xbound(1)=0.0
xbound(2)=0.0
ybound(1)=0.0
ybound(2)=0.0
zbound(1)=0.0
zbound(2)=0.0
delx(1)=0.0
delx(2)=0.0
delx(3)=0.0
lattype='none'
read(5,latcard)
if (lattype .eq. 'none') goto 2
call latgen1()
call defvalue()
goto 1
2 continue
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine latgen1()
include 'createAtoms.h'
character*8 lattype
common /latinfo/ lattype,alat(3),xrot(3),yrot(3),zrot(3),
* periodicity(3),xbound(2),ybound(2),zbound(2),strain(3),
* delx(3),rcell(3),ccell(nelmax)
common /tmpvar/ avec(3),bvec(3),cvec(3),xold(3),yold(3),zold(3),
* avecp(3),bvecp(3),cvecp(3),rcellp(3),roter(3,3),rvs(3,100)
data xold/1.0,0.0,0.0/,yold/0.0,1.0,0.0/,zold/0.0,0.0,1.0/
namelist /subcard/ rcell,ccell
call storelat(lattype,avec,bvec,cvec)
if (xbound(1) .eq. xbound(2)) then
xbound(1)=perlb(1)
xbound(2)=perub(1)
endif
if (ybound(1) .eq. ybound(2)) then
ybound(1)=perlb(2)
ybound(2)=perub(2)
endif
if (zbound(1) .eq. zbound(2)) then
zbound(1)=perlb(3)
zbound(2)=perub(3)
endif
vol=alat(1)*alat(2)*alat(3)
xnorm=sqrt((xrot(1)*alat(2)*alat(3))**2+
* (xrot(2)*alat(1)*alat(3))**2+
* (xrot(3)*alat(1)*alat(2))**2)
ynorm=sqrt((yrot(1)*alat(2)*alat(3))**2+
* (yrot(2)*alat(1)*alat(3))**2+
* (yrot(3)*alat(1)*alat(2))**2)
znorm=sqrt((zrot(1)*alat(2)*alat(3))**2+
* (zrot(2)*alat(1)*alat(3))**2+
* (zrot(3)*alat(1)*alat(2))**2)
periodicity(1)=vol*periodicity(1)/xnorm
periodicity(2)=vol*periodicity(2)/ynorm
periodicity(3)=vol*periodicity(3)/znorm
do 5 i=1,3
xrot(i)=vol*xrot(i)/(alat(i)*xnorm)
yrot(i)=vol*yrot(i)/(alat(i)*ynorm)
zrot(i)=vol*zrot(i)/(alat(i)*znorm)
5 continue
do 10 i=1,3
do 10 j=1,3
roter(i,j)=xold(i)*xrot(j)+yold(i)*yrot(j)+zold(i)*zrot(j)
10 continue
do 20 i=1,3
avecp(i)=0.0
bvecp(i)=0.0
cvecp(i)=0.0
do 20 j=1,3
avecp(i)=avecp(i)+roter(i,j)*avec(j)*0.5*alat(j)
bvecp(i)=bvecp(i)+roter(i,j)*bvec(j)*0.5*alat(j)
cvecp(i)=cvecp(i)+roter(i,j)*cvec(j)*0.5*alat(j)
20 continue
do 30 i=1,3
avec(i)=avecp(i)
bvec(i)=bvecp(i)
cvec(i)=cvecp(i)
30 continue
call applystrain()
nrange=20
nsmall=0
do 300 ii=-nrange,nrange
do 300 jj=-nrange,nrange
do 300 kk=-nrange,nrange
xt=ii*avec(1)+jj*bvec(1)+kk*cvec(1)
if (xt .ge. periodicity(1)-0.1 .or. xt .lt. -0.1) goto 300
yt=ii*avec(2)+jj*bvec(2)+kk*cvec(2)
if (yt .ge. periodicity(2)-0.1 .or. yt .lt. -0.1) goto 300
zt=ii*avec(3)+jj*bvec(3)+kk*cvec(3)
if (zt .ge. periodicity(3)-0.1 .or. zt .lt. -0.1) goto 300
nsmall=nsmall+1
rvs(1,nsmall)=xt
rvs(2,nsmall)=yt
rvs(3,nsmall)=zt
300 continue
nxa=xbound(1)/periodicity(1)-2
nxb=xbound(2)/periodicity(1)+2
nya=ybound(1)/periodicity(2)-2
nyb=ybound(2)/periodicity(2)+2
nza=zbound(1)/periodicity(3)-2
nzb=zbound(2)/periodicity(3)+2
1 rcell(1)=-100.0
read(5,subcard)
if (rcell(1) .eq. -100.0) then
return
endif
do 190 m = 2,ntypes
ccell(m)=min(1.0,ccell(m-1)+ccell(m))
190 continue
do 800 i=1,3
rcellp(i)=0.0
do 800 j=1,3
rcellp(i)=rcellp(i)+roter(i,j)*rcell(j)*alat(j)
800 continue
do 801 i=1,3
rcell(i)=rcellp(i)*(1.0+strain(i))
801 continue
do 250 nn=1,nsmall
do 251 ii=nxa,nxb
xt=ii*periodicity(1)+rvs(1,nn)+rcell(1)
if (xt .ge. xbound(2) .or. xt .lt. xbound(1)) goto 251
do 252 jj=nya,nyb
yt=jj*periodicity(2)+rvs(2,nn)+rcell(2)
if (yt .ge. ybound(2) .or. yt .lt. ybound(1)) goto 252
do 253 kk=nza,nzb
zt=kk*periodicity(3)+rvs(3,nn)+rcell(3)
if (zt .ge. zbound(2) .or. zt .lt. zbound(1)) goto 253
tst=ranl()
it=0
do 201 m=ntypes,1,-1
if (tst .le. ccell(m)) it=m
201 continue
natoms=natoms+1
if (natoms .gt. natmax) then
write(6,9601)natmax
9601 format('number of atoms exceeds maximum dimension: ',i10)
endif
itype(natoms)=it
nhittag(natoms)=0
nntype(it)=nntype(it)+1
rv(1,natoms)=xt+delx(1)
rv(2,natoms)=yt+delx(2)
rv(3,natoms)=zt+delx(3)
if (rv(1,natoms) .lt. perlb(1))
* rv(1,natoms)=rv(1,natoms)+perlen(1)
if (rv(1,natoms) .ge. perub(1))
* rv(1,natoms)=rv(1,natoms)-perlen(1)
if (rv(2,natoms) .lt. perlb(2))
* rv(2,natoms)=rv(2,natoms)+perlen(2)
if (rv(2,natoms) .ge. perub(2))
* rv(2,natoms)=rv(2,natoms)-perlen(2)
if (rv(3,natoms) .lt. perlb(3))
* rv(3,natoms)=rv(3,natoms)+perlen(3)
if (rv(3,natoms) .ge. perub(3))
* rv(3,natoms)=rv(3,natoms)-perlen(3)
253 continue
252 continue
251 continue
250 continue
goto 1
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine defvalue()
include 'createAtoms.h'
integer oldtype
character*8 lattype
common /latinfo/ lattype,alat(3),xrot(3),yrot(3),zrot(3),
* periodicity(3),xbound(2),ybound(2),zbound(2),strain(3),
* delx(3),rcell(3),ccell(nelmax)
common /tmpvar/ avec(3),bvec(3),cvec(3),xold(3),yold(3),zold(3),
* avecp(3),bvecp(3),cvecp(3),rcellp(3),roter(3,3),rvs(3,100)
dimension cent(3),plane(3),planec(3),surface(3)
namelist /defcard/ xmin,xmax,ymin,ymax,zmin,zmax,
* oldtype,newtype,prob,cent,plane,dist,surface,
* nramp,ndirec,dismax,theta,phi,nscrew,anglemin,anglemax,
* xbottom,xheight,zbottom,zheight,nrepeatsx,nrepeatsz
data pi/3.1415927/
1 continue
xmin=-1.0e-20
xmax=1.0e20
ymin=-1.0e20
ymax=1.0e20
zmin=-1.0e20
zmax=1.0e20
oldtype=0
dist=-1.0
newtype=-11
cent(1)=0.5*(perlb(1)+perub(1))
cent(2)=0.5*(perlb(2)+perub(2))
cent(3)=0.5*(perlb(3)+perub(3))
surface(1)=0.0
surface(2)=0.0
surface(3)=0.0
nramp=1
nscrew=0
anglemin=0.0
anglemax=360.0
ndirec=3
dismax=0.0
theta=5000.0
phi=5000.0
nrepeatsx=0
nrepeatsz=0
read(5,defcard)
if (newtype .lt. -10 .and. dismax .eq. 0.0) return
if (newtype .gt. ntypes) then
ntypes=ntypes+1
if (oldtype .ne. 0) then
amass(ntypes)=amass(oldtype)
else
amass(ntypes)=amass(ntypes-1)
endif
endif
if (surface(1) .eq. 1.0) then
perlb(1)=perlb(1)-10.0
perub(1)=perub(1)+10.0
perlen(1)=perub(1)-perlb(1)
endif
if (surface(2) .eq. 1.0) then
perlb(2)=perlb(2)-10.0
perub(2)=perub(2)+10.0
perlen(2)=perub(2)-perlb(2)
endif
if (surface(3) .eq. 1.0) then
perlb(3)=perlb(3)-10.0
perub(3)=perub(3)+10.0
perlen(3)=perub(3)-perlb(3)
endif
vol=alat(1)*alat(2)*alat(3)
if (dismax .ne. 0.0) then
astart=anglemin
afinish=anglemax
if (anglemax .lt. anglemin) then
anglemin=afinish
anglemax=astart
endif
if (nramp .eq. 1) then
pstart=xmin
pfinish=xmax
if (xmax .lt. xmin) then
xmin=pfinish
xmax=pstart
endif
endif
if (nramp .eq. 2) then
pstart=ymin
pfinish=ymax
if (ymax .lt. ymin) then
ymin=pfinish
ymax=pstart
endif
endif
if (nramp .eq. 3) then
pstart=zmin
pfinish=zmax
if (zmax .lt. zmin) then
zmin=pfinish
zmax=pstart
endif
endif
do 33 i=1,natoms
if (rv(1,i) .lt. xmin) goto 33
if (rv(1,i) .gt. xmax) goto 33
if (rv(2,i) .lt. ymin) goto 33
if (rv(2,i) .gt. ymax) goto 33
if (rv(3,i) .lt. zmin) goto 33
if (rv(3,i) .gt. zmax) goto 33
if (nscrew .eq. 0) then
if (nramp .eq. 0) then
rv(ndirec,i)=rv(ndirec,i)+dismax
else
rv(ndirec,i)=rv(ndirec,i)+
* (rv(nramp,i)-pstart)/(pfinish-pstart)*dismax
endif
else
if (ndirec .eq. 1) then
dx=rv(2,i)-cent(2)
dy=rv(3,i)-cent(3)
dxmax=abs(ymax-cent(2))
dxmin=abs(ymin-cent(2))
rmax=min(dxmin,dxmax)
endif
if (ndirec .eq. 2) then
dx=rv(3,i)-cent(3)
dy=rv(1,i)-cent(1)
dxmax=abs(zmax-cent(3))
dxmin=abs(zmin-cent(3))
rmax=min(dxmin,dxmax)
endif
if (ndirec .eq. 3) then
dx=rv(1,i)-cent(1)
dy=rv(2,i)-cent(2)
dxmax=abs(xmax-cent(1))
dxmin=abs(xmin-cent(1))
rmax=min(dxmin,dxmax)
endif
dr=sqrt(dx**2+dy**2)
if (dr .eq. 0.0) then
angle=0.0
else
angle=acos(dx/dr)*180.0/pi
if (dy .lt. 0.0) angle=360.0-angle
endif
if (angle .lt. anglemin) goto 33
if (angle .gt. anglemax) goto 33
if (nscrew .eq. 1) then
astart1=astart
afinish1=afinish
else if (nscrew .eq. 2) then
if (astart .eq. 0.0 .and.
* afinish .eq. 180.0) then
if (dr .le. dxmax) then
astart1=astart
else
astart1=acos(dxmax/dr)*180.0/pi
endif
if (dr .le. dxmin) then
afinish1=afinish
else
afinish1=180.0-acos(dxmin/dr)*180.0/pi
endif
else if (astart .eq. 180.0 .and.
* afinish .eq. 0.0) then
if (dr .le. dxmin) then
astart1=astart
else
astart1=180.0-acos(dxmin/dr)*180.0/pi
endif
if (dr .le. dxmax) then
afinish1=afinish
else
afinish1=acos(dxmax/dr)*180.0/pi
endif
else if (astart .eq. 180.0 .and.
* afinish .eq. 360.0) then
if (dr .le. dxmin) then
astart1=astart
else
astart1=180.0+acos(dxmin/dr)*180.0/pi
endif
if (dr .le. dxmax) then
afinish1=afinish
else
afinish1=360.0-acos(dxmax/dr)*180.0/pi
endif
else if (astart .eq. 360.0 .and.
* afinish .eq. 180.0) then
if (dr .le. dxmax) then
astart1=astart
else
astart1=360.0-acos(dxmax/dr)*180.0/pi
endif
if (dr .le. dxmin) then
afinish1=afinish
else
afinish1=180.0+acos(dxmin/dr)*180.0/pi
endif
else
write(6,*)'problem in nscrew=2'
stop
endif
endif
if (dr .gt. rmax) dr=rmax
rv(ndirec,i)=rv(ndirec,i)+dr/rmax*
* (angle-astart1)/(afinish1-astart1)*dismax
endif
33 continue
else if (theta .eq. 5000.0 .and. phi .eq. 5000.0 .and.
* nrepeatsx*nrepeatsz .eq. 0.0) then
if (dist .lt. 0.0) then
do 2 i=1,natoms
if (rv(1,i) .lt. xmin) goto 2
if (rv(1,i) .gt. xmax) goto 2
if (rv(2,i) .lt. ymin) goto 2
if (rv(2,i) .gt. ymax) goto 2
if (rv(3,i) .lt. zmin) goto 2
if (rv(3,i) .gt. zmax) goto 2
if (oldtype .eq. itype(i) .or. oldtype .eq. 0) then
if (ranl() .lt. prob) itype(i)=newtype
endif
2 continue
else
if (theta .eq. 5000.0) then
pnorm=sqrt((plane(1)*alat(2)*alat(3))**2+
* (plane(2)*alat(1)*alat(3))**2+
* (plane(3)*alat(1)*alat(2))**2)
do 30 n=1,3
plane(n)=vol*plane(n)/(alat(n)*pnorm)
30 continue
do 20 n=1,3
planec(n)=0.0
do 20 m=1,3
planec(n)=planec(n)+roter(n,m)*plane(m)
20 continue
else
planec(1)=sin(pi/180.0*theta)*cos(pi/180.0*phi)
planec(3)=sin(pi/180.0*theta)*sin(pi/180.0*phi)
planec(2)=cos(pi/180.0*theta)
endif
do 22 i=1,natoms
dis1=rv(1,i)-cent(1)
dis2=rv(2,i)-cent(2)
dis3=rv(3,i)-cent(3)
if (dis1*planec(1)+dis2*planec(2)+dis3*planec(3)
* .gt. dist) then
if (oldtype .eq. itype(i) .or. oldtype .eq. 0)
* itype(i)=newtype
endif
22 continue
endif
else if (nrepeatsx*nrepeatsz .ne. 0) then
do 4 i=1,natoms
if (rv(1,i) .le. xmin .or. rv(1,i) .ge. xmax) then
ylocalx = xbottom
else
ylocalx = xbottom + 0.5*xheight*(1.0+
* cos(2.0*pi*nrepeatsx/(xmax-xmin)*
* (rv(1,i)-(xmin+xmax)/(2.0*nrepeatsx))))
endif
if (rv(3,i) .le. zmin .or. rv(3,i) .ge. zmax) then
ylocalz = zbottom
else
ylocalz = zbottom + 0.5*zheight*(1.0+
* cos(2.0*pi*nrepeatsz/(zmax-zmin)*
* (rv(3,i)-(zmin+zmax)/(2.0*nrepeatsz))))
endif
ylocal = ylocalx
if (ylocal .gt. ylocalz) ylocal = ylocalz
if (rv(2,i) .gt. ylocal) then
itype(i)=newtype
endif
4 continue
endif
goto 1
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine storelat(latty,avec,bvec,cvec)
character*8 lstored(10),latty
dimension avecs(3,10),bvecs(3,10),cvecs(3,10)
dimension avec(3),bvec(3),cvec(3)
data lstored/'fcc ','bcc ','sc ',7*' '/
data avecs/1.0,1.0,0.0, 1.0,1.0,-1.0, 2.0,0.0,0.0, 21*0.0/
data bvecs/0.0,1.0,1.0, -1.0,1.0,1.0, 0.0,2.0,0.0, 21*0.0/
data cvecs/1.0,0.0,1.0, 1.0,-1.0,1.0, 0.0,0.0,2.0, 21*0.0/
imatch=0
do 10 i=1,10
if (latty .ne. lstored(i)) goto 10
imatch = 1
do 5 j=1,3
avec(j)=avecs(j,i)
bvec(j)=bvecs(j,i)
cvec(j)=cvecs(j,i)
5 continue
10 continue
if (imatch .eq. 1) return
write(6,15) latty
15 format('could not find lattice type ')
stop
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine applystrain()
include 'createAtoms.h'
character*8 lattype
common /latinfo/ lattype,alat(3),xrot(3),yrot(3),zrot(3),
* periodicity(3),xbound(2),ybound(2),zbound(2),strain(3),
* delx(3),rcell(3),ccell(nelmax)
common /tmpvar/ avec(3),bvec(3),cvec(3),xold(3),yold(3),zold(3),
* avecp(3),bvecp(3),cvecp(3),rcellp(3),roter(3,3),rvs(3,100)
do i=1,3
periodicity(i)=(1.0+strain(i))*periodicity(i)
avec(i)=(1.0+strain(i))*avec(i)
bvec(i)=(1.0+strain(i))*bvec(i)
cvec(i)=(1.0+strain(i))*cvec(i)
enddo
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine hitvalue()
include 'createAtoms.h'
namelist /hitcard/ xmin,xmax,ymin,ymax,zmin,zmax,nhit
nhitcards=0
1001 nhit=0
xmin=-10000.0
xmax=10000.0
ymin=-10000.0
ymax=10000.0
zmin=-10000.0
zmax=10000.0
read(5,hitcard)
if (nhit .eq. 0) then
return
endif
nchange=0
do 1002 i=1,natoms
if (rv(1,i) .lt. xmin) goto 1002
if (rv(1,i) .gt. xmax) goto 1002
if (rv(2,i) .lt. ymin) goto 1002
if (rv(2,i) .gt. ymax) goto 1002
if (rv(3,i) .lt. zmin) goto 1002
if (rv(3,i) .gt. zmax) goto 1002
nchange=nchange+1
nhittag(i)=nhit
1002 continue
if (nchange .ne. 0) nhitcards=nhitcards+1
goto 1001
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine systemshift()
include 'createAtoms.h'
namelist /shiftcard/ mode
mode=0
read(5,shiftcard)
if (mode .eq. 0) return
if (mode .eq. 1) then
top=perub(2)
bottom=perlb(2)
sdely=0.0
do i=1,natoms
if (nhittag(i) .eq. 3 .and. top .gt. rv(2,i)) top=rv(2,i)
if (nhittag(i) .eq. 4 .and. bottom .lt. rv(2,i))
* bottom=rv(2,i)
enddo
sdely=-0.5*(top+bottom)
do i=1,natoms
rv(2,i)=rv(2,i)+sdely
enddo
endif
if (mode .eq. 2) then
c This shift only assumes a maximum of two plane spacings.
xmin=10.0
do i=1,natoms
if (rv(1,i) .lt. xmin) xmin=rv(1,i)
enddo
shiftneg=-10.0
shiftpos=10.0
do i=1,natoms
shift=rv(1,i)-xmin
shift=shift-perlen(1)*nint(shift/perlen(1))
if (shift .gt. 0.01 .and. shift .lt. shiftpos)
* shiftpos=shift
if (shift .lt. -0.01 .and. shift .gt. shiftneg)
* shiftneg=shift
enddo
if (-shiftneg .gt. shiftpos) then
shift=0.5*shiftneg
else
shift=0.5*shiftpos
endif
perlb(1)=xmin+shift
perub(1)=perlb(1)+perlen(1)
endif
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine colect()
include 'createAtoms.h'
do i=1,natoms
if (rv(1,i) .lt. perlb(1)) rv(1,i)=rv(1,i)+perlen(1)
if (rv(1,i) .ge. perub(1)) rv(1,i)=rv(1,i)-perlen(1)
if (rv(2,i) .lt. perlb(2)) rv(2,i)=rv(2,i)+perlen(2)
if (rv(2,i) .ge. perub(2)) rv(2,i)=rv(2,i)-perlen(2)
if (rv(3,i) .lt. perlb(3)) rv(3,i)=rv(3,i)+perlen(3)
if (rv(3,i) .ge. perub(3)) rv(3,i)=rv(3,i)-perlen(3)
enddo
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
real function ranl()
common /iran/ iseed
iseed=iseed*125
iseed=iseed-2796203*(iseed/2796203)
ranl=abs(iseed/2796203.0)
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine delete()
include 'createAtoms.h'
nw_del=0
do i=1,ntypes
nntype(i)=0
enddo
natoms0=natoms
natoms=0
do 1 i=1,natoms0
if (itype(i) .eq. -1) nw_del=nw_del+1
if (itype(i) .le. 0) goto 1
natoms=natoms+1
ntag(natoms)=natoms
itype(natoms)=itype(i)
nhittag(natoms)=nhittag(i)
nntype(itype(natoms))=nntype(itype(natoms))+1
rv(1,natoms)=rv(1,i)
rv(2,natoms)=rv(2,i)
rv(3,natoms)=rv(3,i)
1 continue
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine randomize()
include 'createAtoms.h'
dimension rv0(6,natmax),itype0(natmax),nhittag0(natmax)
do 1 i=1,natoms
itype0(i)=itype(i)
if (nhitcards .gt. 0) nhittag0(i)=nhittag(i)
rv0(1,i)=rv(1,i)
rv0(2,i)=rv(2,i)
rv0(3,i)=rv(3,i)
rv0(4,i)=rv(4,i)
rv0(5,i)=rv(5,i)
rv0(6,i)=rv(6,i)
1 continue
im=1
10 im=im*2
if (im .lt. natoms) goto 10
ia=365
ic=8161
ivalue=mod(ilatseed,natoms)
do 2 i=1,natoms
20 ivalue=mod(ia*ivalue+ic,im)
if (ivalue .ge. natoms) goto 20
ntmp=ivalue+1
if (ntmp .ge. 1 .and. ntmp .le. natoms) then
ntag(ntmp)=i
itype(ntmp)=itype0(i)
if (nhitcards .gt. 0) nhittag(ntmp)=nhittag0(i)
rv(1,ntmp)=rv0(1,i)
rv(2,ntmp)=rv0(2,i)
rv(3,ntmp)=rv0(3,i)
rv(4,ntmp)=rv0(4,i)
rv(5,ntmp)=rv0(5,i)
rv(6,ntmp)=rv0(6,i)
endif
2 continue
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine write0(dumpfile)
include 'createAtoms.h'
character*32 dumpfile, typ
open (unit=2,file=dumpfile)
write(2,*) natoms
write(2,100) 'GaN lattice'
do i=1,natoms
if (itype(i).eq.1) then
typ = 'Ga'
else
typ = 'N'
endif
write(2,150) typ(1:length(typ)),rv(1,i),rv(2,i),rv(3,i)
enddo
50 format(i8)
100 format(a)
150 format(a3,3f13.5)
close(2)
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine write1(dumpfile)
include 'createAtoms.h'
character*32 dumpfile
zero=0.0
open (unit=2,file=dumpfile)
write(2,*) 'ITEM: NUMBER OF ATOMS'
write(2,*) natoms
write(2,*) 'ITEM: BOX BOUNDS'
write(2,*) perlb(1),perub(1)
write(2,*) perlb(2),perub(2)
write(2,*) perlb(3),perub(3)
write(2,*) 'ITEM: TIME'
write(2,*) 0.0
write(2,*) 'ITEM: ATOMS'
if (nhitcards .eq. 0) then
do i=1,natoms
write(2,*) ntag(i),itype(i),rv(1,i),rv(2,i),rv(3,i)
enddo
else
do i=1,natoms
write(2,*) ntag(i),itype(i),nhittag(i),
* rv(1,i),rv(2,i),rv(3,i)
enddo
endif
close(2)
open(unit=2,file=dumpfile(1:index(dumpfile,' ')-1)//'.vel')
write(2,*) 'ITEM: NUMBER OF ATOMS'
write(2,*) natoms
write(2,*) 'ITEM: BOUNDARY VELOCITY'
write(2,*) zero,zero,zero
write(2,*) 'ITEM: TIME'
write(2,*) zero
write(2,*) 'ITEM: VELOCITIES'
do i=1,natoms
write(2,*) ntag(i),rv(4,i),rv(5,i),rv(6,i)
enddo
close(2)
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine write2(dumpfile)
include 'createAtoms.h'
data conmas/1.0365e-4/
character*32 dumpfile
open (unit=2,file=dumpfile)
zero=0.0
write(2,*) 'using dynamo'
write(2,9502) natoms,ntypes,zero
9502 format(2i10,e15.8)
write(2,9503) (perub(i),i=1,3),(perlb(i),i=1,3)
9503 format(3e25.16)
write(2,9504) (amass(i)*conmas,ielement(i),i=1,ntypes)
9504 format(e25.16,i10)
write(2,9505) ((rv(i,j),i=1,6),itype(j),j=1,natoms)
9505 format(3e25.16/3e25.16/i10)
write(2,9503) zero,zero,zero
close(2)
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine write3(dumpfile)
include 'createAtoms.h'
character*32 dumpfile
open (unit=2,file=dumpfile)
area=(1.0-float(nw_del)/float(natoms0))*perlen(2)*perlen(3)
write(2,*) 'Start File for LAMMPS ',area
write(2,*)
write(2,*) natoms,' atoms'
write(2,*)
write(2,*) ntypes,' atom types'
write(2,*)
write(2,*) perlb(1),perub(1),' xlo xhi'
write(2,*) perlb(2),perub(2),' ylo yhi'
write(2,*) perlb(3),perub(3),' zlo zhi'
if (xy .ge. 1.0e8 .or. xz .ge. 1.0e8 .or. yz .ge. 1.0e8) then
write(2,*)
else
write(2,*)
write(2,*)xy,xz,yz,' xy xz yz'
write(2,*)
endif
write(2,*) 'Masses'
write(2,*)
do i=1,ntypes
write(2,*)i,amass(i)
enddo
write(2,*)
write(2,*) 'Atoms'
write(2,*)
do i=1,natoms
write(2,*) i,itype(i),rv(1,i),rv(2,i),rv(3,i)
enddo
write(2,*)
write(2,*) 'Velocities'
write(2,*)
zero=0.0
do i=1,natoms
write(2,*) i,rv(4,i),rv(5,i),rv(6,i)
enddo
close(2)
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine write4(dumpfile)
include 'createAtoms.h'
character*32 dumpfile,full
dumpfile='ens'
full=dumpfile(1:index(dumpfile,' ')-1)//'.case'
open(unit=2,file=full)
write(2,101)
101 format('FORMAT')
write(2,102)
102 format('type: ensight')
write(2,103)
103 format('GEOMETRY')
write(2,104)'ens.case.****.geo'
104 format('model: 1 ',a30)
write(2,105)
105 format('VARIABLE')
write(2,111)
111 format('TIME')
write(2,112)
112 format('time set: 1')
write(2,113)
113 format('number of steps: 1')
write(2,114)
114 format('filename start number: 1')
write(2,115)
115 format('filename increment: 1')
write(2,116)
116 format('time values:')
write(2,117)
117 format('1')
close(2)
full=dumpfile(1:index(dumpfile,' ')-1)//'.case.0001.geo'
open(unit=2,file=full)
write(2,201)
write(2,202)
write(2,203)
write(2,204)
write(2,205)
write(2,206)natoms
do j=1,natoms
write(2,207)j,rv(1,j),rv(2,j),rv(3,j)
enddo
do j=1,ntypes
write(2,208)j
write(2,209)j
write(2,210)
write(2,206)nntype(j)
do k=1,natoms
if (itype(k) .eq. j) write(2,211)k,k
enddo
enddo
201 format('spparks input')
202 format('geometry for element group 1')
203 format('node id given')
204 format('element id given')
205 format('coordinates')
206 format(i8)
207 format(i8,3e12.5)
208 format('part',i2)
209 format('atombox',i2)
210 format('point')
211 format(2i8)
close(2)
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine write5(dumpfile)
include 'createAtoms.h'
character*32 dumpfile,full
character*200 sentence
full=dumpfile(1:index(dumpfile,' ')-1)//'.lat'
open (unit=2,file=full)
call neigen()
nmax=numneigh(1)
do i=1,natoms
if (nmax .lt. numneigh(i)) nmax=numneigh(i)
enddo
write(2,*) 'spparks input lattice'
write(2,*)
write(2,*) '3 dimension'
write(2,*) natoms,' vertices'
write(2,*) nmax,' max connectivity'
write(2,*) perlb(1),perub(1),' xlo xhi'
write(2,*) perlb(2),perub(2),' ylo yhi'
write(2,*) perlb(3),perub(3),' zlo zhi'
write(2,*)
write(2,*) 'Vertices'
write(2,*)
do i=1,natoms
write(2,*) i,rv(1,i),rv(2,i),rv(3,i)
enddo
write(2,*)
write(2,*) 'Edges'
write(2,*)
do i=1,natoms
write(sentence,*) i,(neigh(k,i),k=1,numneigh(i))
write(2,*)sentence(1:len_trim(sentence))
enddo
close(2)
full=dumpfile(1:index(dumpfile,' ')-1)//'.conf'
open (unit=2,file=full)
nvalues=1
write(2,*) 'spparks input configuration'
write(2,*) natoms,nvalues
write(2,*)
do i=1,natoms
write(2,*) i,itype(i)
enddo
close(2)
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine disturb()
include 'createAtoms.h'
dimension strain(3)
namelist /disturbcard/ dismax,strain
dismax=0.0
strain(1)=0.0
strain(2)=0.0
strain(3)=0.0
read(5,disturbcard)
if (dismax .ne. 0.0) then
do 1 i=1,natoms
rv(1,i)=rv(1,i)+dismax*(ranl()-0.5)
rv(2,i)=rv(2,i)+dismax*(ranl()-0.5)
rv(3,i)=rv(3,i)+dismax*(ranl()-0.5)
1 continue
endif
if (strain(1)**2+strain(2)**2+strain(3)**2 .ne. 0.0) then
do 2 i=1,3
perlen(i)=perlen(i)*(1.0+strain(i))
perub(i)=perlb(i)+perlen(i)
2 continue
do 3 i=1,natoms
do 3 j=1,3
rv(j,i)=perlb(j)+(rv(j,i)-perlb(j))*(1.0+strain(j))
3 continue
endif
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine outputfiles()
character*32 dynamo,paradyn,lammps,xyz,ensight,spparks
namelist /filecard/ dynamo,paradyn,lammps,xyz,ensight,spparks
xyz="none"
paradyn="none"
dynamo="none"
lammps="none"
ensight="none"
spparks="none"
read(5,filecard)
if (xyz .ne. "none") call write0(xyz)
if (paradyn .ne. "none") call write1(paradyn)
if (dynamo .ne. "none") call write2(dynamo)
if (lammps .ne. "none") call write3(lammps)
if (ensight .ne. "none") call write4(ensight)
if (spparks .ne. "none") call write5(spparks)
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
integer function length(str)
character*(*) str
n = len(str)
10 if (n.gt.0) then
if (str(n:n).eq.' ') then
n = n - 1
goto 10
endif
endif
length = n
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine velgen()
include 'createAtoms.h'
data boltz/8.617e-5/,conmas/1.0365e-4/
namelist /velcard/ Tbnd,Tmid,naxis,nmesh,equal,ensureTave,
* steeper
dimension Tmesh(200),pmesh(200),vcm(3,200),tmass(200),ekin(200),
* scale(200),ncount(200)
Tbnd=0.0
Tmid=0.0
naxis=1
nmesh=100
equal=1.0
ensureTave=0.0
steeper=0.0
read(5,velcard)
Tbnd=equal*Tbnd
Tmid=equal*Tmid
Tave=0.5*(Tbnd+Tmid)
Tbnd=Tave+(Tbnd-Tave)*(1.0+steeper)
Tmid=Tave+(Tmid-Tave)*(1.0+steeper)
if (Tbnd+Tmid .le. 0.0) then
do i=1,natoms
rv(4,i)=0.0
rv(5,i)=0.0
rv(6,i)=0.0
enddo
else
wmesh=perlen(naxis)/nmesh
do i=1,nmesh
pmesh(i)=(i-0.5)*wmesh
Tmesh(i)=Tmid+(Tbnd-Tmid)*abs(pmesh(i)-0.5*perlen(naxis))/
* (0.5*perlen(naxis))
vcm(1,i)=0.0
vcm(2,i)=0.0
vcm(3,i)=0.0
tmass(i)=0.0
ekin(i)=0.0
ncount(i)=0
enddo
do i=1,natoms
index=(rv(naxis,i)-perlb(naxis))/wmesh+1
it=itype(i)
vnorm=sqrt(Tmesh(index)*boltz/(amass(it)*conmas))
do j=1,3
rv(j+3,i)=vnorm*rgauss()
vcm(j,index)=vcm(j,index)+amass(it)*conmas*rv(j+3,i)
enddo
tmass(index)=tmass(index)+amass(it)*conmas
ncount(index)=ncount(index)+1
enddo
do i=1,nmesh
do j=1,3
vcm(j,i)=vcm(j,i)/tmass(i)
enddo
enddo
do i=1,natoms
index=(rv(naxis,i)-perlb(naxis))/wmesh+1
do j=1,3
rv(j+3,i)=rv(j+3,i)-vcm(j,index)
enddo
enddo
do i=1,natoms
index=(rv(naxis,i)-perlb(naxis))/wmesh+1
it=itype(i)
ekin(index)=ekin(index)+0.5*amass(it)*conmas*
* (rv(4,i)**2+rv(5,i)**2+rv(6,i)**2)
enddo
do i=1,nmesh
treal=2.0/(3.0*boltz)*ekin(i)/ncount(i)
if (treal .eq. 0.0) then
scale(i)=0.0
else
scale(i)=sqrt(Tmesh(i)/treal)
endif
enddo
do i=1,natoms
index=(rv(naxis,i)-perlb(naxis))/wmesh+1
do j=4,6
rv(j,i)=rv(j,i)*scale(index)
enddo
enddo
if (ensureTave .eq. 1.0) then
ekint=0.0
do i=1,natoms
it=itype(i)
ekint=ekint+0.5*amass(it)*conmas*
* (rv(4,i)**2+rv(5,i)**2+rv(6,i)**2)
enddo
trealt=2.0/(3.0*boltz)*ekint/natoms
if (trealt .eq. 0.0) then
scalet=0.0
else
scalet=sqrt(Tave/trealt)
endif
do i=1,natoms
do j=4,6
rv(j,i)=rv(j,i)*scalet
enddo
enddo
endif
endif
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
function rgauss()
data twopi/6.283185308/
1 continue
a1 = ranl()
if (a1 .eq. 0.0) goto 1
a2 = ranl()
rgauss = sqrt(-2.0*log(a1))*cos(twopi*a2)
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine neigen()
include 'createAtoms.h'
dimension idbox(100,100,100,100),nbox(100,100,100)
cutoff=100.0
xmin=rv(1,1)
xmax=rv(1,1)
ymin=rv(2,1)
ymax=rv(2,1)
zmin=rv(3,1)
zmax=rv(3,1)
do i=2,natoms
if (xmin .gt. rv(1,i)) xmin=rv(1,i)
if (xmax .lt. rv(1,i)) xmax=rv(1,i)
if (ymin .gt. rv(2,i)) ymin=rv(2,i)
if (ymax .lt. rv(2,i)) ymax=rv(2,i)
if (zmin .gt. rv(3,i)) zmin=rv(3,i)
if (zmax .lt. rv(3,i)) zmax=rv(3,i)
dx=rv(1,i)-rv(1,1)
dy=rv(2,i)-rv(2,1)
dz=rv(3,i)-rv(3,1)
dx=dx-perlen(1)*nint(dx/perlen(1))
dy=dy-perlen(2)*nint(dy/perlen(2))
dz=dz-perlen(3)*nint(dz/perlen(3))
tmp=dx**2+dy**2+dz**2
if (tmp .lt. cutoff) cutoff=tmp
enddo
cutoff=sqrt(cutoff)+0.1
nx=(xmax-xmin)/cutoff+1
ny=(ymax-ymin)/cutoff+1
nz=(zmax-zmin)/cutoff+1
if (nx .gt. 100) nx=100
if (ny .gt. 100) ny=100
if (nz .gt. 100) nz=100
delx=(xmax-xmin)/nx
dely=(ymax-ymin)/ny
delz=(zmax-zmin)/nz
do n1=1,nx
do n2=1,ny
do n3=1,nz
nbox(n1,n2,n3)=0
enddo
enddo
enddo
do j=1,natoms
numneigh(j)=0
n1=(rv(1,j)-xmin)/delx+1
if (n1 .lt. 1) n1=1
if (n1 .gt. nx) n1=nx
n2=(rv(2,j)-ymin)/dely+1
if (n2 .lt. 1) n2=1
if (n2 .gt. ny) n2=ny
n3=(rv(3,j)-zmin)/delz+1
if (n3 .lt. 1) n3=1
if (n3 .gt. nz) n3=nz
nbox(n1,n2,n3)=nbox(n1,n2,n3)+1
idbox(n1,n2,n3,nbox(n1,n2,n3))=j
enddo
do j=1,natoms
n1=(rv(1,j)-xmin)/delx+1
if (n1 .lt. 1) n1=1
if (n1 .gt. nx) n1=nx
n2=(rv(2,j)-ymin)/dely+1
if (n2 .lt. 1) n2=1
if (n2 .gt. ny) n2=ny
n3=(rv(3,j)-zmin)/delz+1
if (n3 .lt. 1) n3=1
if (n3 .gt. nz) n3=nz
do na=n1-1,n1+1
if (na .lt. 1) then
naa=nx
else if (na .gt. nx) then
naa=1
else
naa=na
endif
do nb=n2-1,n2+1
if (nb .lt. 1) then
nbb=ny
else if (nb .gt. ny) then
nbb=1
else
nbb=nb
endif
do nc=n3-1,n3+1
if (nc .lt. 1) then
ncc=nz
else if (nc .gt. nz) then
ncc=1
else
ncc=nc
endif
do 100 ndd=1,nbox(naa,nbb,ncc)
k=idbox(naa,nbb,ncc,ndd)
if (k .ge. j) goto 100
dx=rv(1,k)-rv(1,j)
dy=rv(2,k)-rv(2,j)
dz=rv(3,k)-rv(3,j)
dx=dx-perlen(1)*nint(dx/perlen(1))
dy=dy-perlen(2)*nint(dy/perlen(2))
dz=dz-perlen(3)*nint(dz/perlen(3))
tmp=sqrt(dx**2+dy**2+dz**2)
if (tmp .le. cutoff) then
numneigh(j)=numneigh(j)+1
neigh(numneigh(j),j)=k
numneigh(k)=numneigh(k)+1
neigh(numneigh(k),k)=j
endif
100 continue
enddo
enddo
enddo
enddo
return
end
c4science · Help