Skip to content

GitLab

P
PolyBuild
Commit cde68ad2 authored by dundaryilmaz's avatar dundaryilmaz
Browse files
Options

Initial commit

parents
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 1777 additions and 0 deletions
Docs/buildingamorphouschains.txt 0 → 100644
View file @ cde68ad2
Monomer File
Similar to unitcells of crystal structures, polymers composed of monomers, thus we need to supply monomer information to polybuild.
Monomers are building blocks of polymers. We can think of beads as monomers.
Monomer file
Docs/howtobuildcrystalstructures.txt 0 → 100644
View file @ cde68ad2
Building Crystal Structures
Building crystal structures fairly straigth forward. There are two steps:
First step you need to prepare a unitcell file which contains definition of conventional unit cell of the desired material.
Format of this file is similar to VASP or ABINIT's POSCAR file:
12
1.0
6.46600 0.000000 0.000
0.00000 5.599720 0.000
0.00000 0.000000 8.075
Direct
3 0.125000 0.850000 0.300000 3 S 59
3 0.375000 0.350000 0.300000 3 S 59
3 0.625000 0.850000 0.300000 3 S 59
3 0.875000 0.350000 0.300000 3 S 59
...
...
Firs line defines the number of atoms in the unit cell. Second line is a scaling factor which modifies lattice constant. Default value is 1.0. Following three lines define lattice constants (i.e. h-matrix) of the unit cell. The unit is in Angstroms.
Next line defines the format for the positions of atoms in the unit cell. Direct corresponds to fractional coordinates. Currently only fractional coordinates are supported. The conversion between fractional and cartesian is:
R_c= h * R_f
Here h stands for h-matrix and R_f and R_c are fractional and cartesian coordinates.
Following lines define positions of atoms. First column is the atom type. Columns 2 to 5 are fractional coordinates. 6th column is again atom type. 7th and 8th columns are element symbol and mass respectively. Here atom type definition is for LAMMPS only. LAMMPS labels atoms with type indices defined at the input script. For example if you are running a simulation which contains Mo , O, S atoms you may label these three elements as 1,2,3 in the force field definition. Then you need to set same labels in your unit cell file in order to generate compatible format.
Defining Crystal Regions in Simulation Box
Defining regions to be filled with crystal structures within the simulation box set up with "blocks" keyword at the input script.
blocks 2
0.0 0.0 10.0 52.0 64.0 28.0 2.0 mos2.ucell
0.0 0.0 50.0 52.0 64.0 15.0 2.0 graphene.ucell
Here number following blocks keyword defines number of regions to be filled with crystal structures. The definitions of these regions followed by this line. With this setup after reading blocks keyword, polybuild will expect two lines to setup crystal regions. Currently polybuild support only rectangular prisms as crystal regions. First three floating numbers define the lower left corner of the region. Next three numbers define length of the region in each direction. Finally last floating point number defines width of the buffer region between the crystal region and polymer structures. According to this example polybuild will avoid placing polymers within 2 A buffer region as well as the crystal structure regions. Buffer width followed by the file name which contains unitcell definitions.
Lattice Mismatch
Starting from the origin of each box, polybuild places unitcells defined in the unitcell file, to fill the region. Number of unit cells in each direction calculated by dividing box length to lattice constant:
N_unitcell_x=NINT(Boxlength_x/h_xx)
If boxlength is integer multiple of lattice constant polybuild will place accordingly. However if it is not polybuild will match the lattice constant to fit into box. For example lets assume box length in x direction is 110A and the lattice constant h_xx=9.8A.
NINT ( nearest integer ) function will yield:
N_unitcell_x=NINT(110/9.8)=NINT(11.22)=11
Polybuild will fit 11 unitcells in x direction by changing lattice constant accordingly:
h_xx_new=Boxlength_x/N_unitcell_x=110/11=10A.
This corresponds to a tensile strain applied in x direction with amount of e_x=(10.0-9.8)/9.8*100= %2.
This feature is really helpfull in building multiple crystal structures as layered configuration. Assume that you are building a composite material consists of two layers of crystal regions. Layers stacked in z direction and periodic boundary conditions applied in x and y directions. In this case you need to find common dimensions in x and y directions to have minumum amount of strain on each layer. For such systems polybuild automatically adjust unitcells by calculating required stress on each structure.
Docs/howtocompile.txt 0 → 100644
View file @ cde68ad2
Polybuild User Manual
1- Compile and Install
Compiling chrpak and geometry libraries.
There are two external libraries requierd to compile polybuild. source code is also located at the Libraries folder.
If lcharpak.a and lgeometry.a files are not compatible with your system you need to compile these two libraries from the source code.
Creating lchrpak.a and lgeometry.a:
To create the library file first you need to compile f90split.f90
gfortran f90split.f90 -o f90split
Make the shell script chrpak.sh executable:
chmod +x chrpak.sh
Run the shell script :
./chrpak.sh
Creating lgeometry.a is similar. You need to run shell script geometry.sh this time.
chmod +x geometry.sh
./geometry.sh
Compiling polybuild is easy. Default compiler is Intel Fortran. In the makefile setting
F90=gfortran
will set the GNU/Gfortran as the Fortran compiler. For best results gcc version > 5.0 recommended.
Change the directory to Src
cd Src
Compile the source codes with:
make.
If everything goes well you should have the executable polybuild.
LICENSE 0 → 100644
View file @ cde68ad2
This diff is collapsed.
Libraries/chrpak.f90 0 → 100644
View file @ cde68ad2
This diff is collapsed.
Libraries/chrpak.sh 0 → 100755
View file @ cde68ad2
#!/bin/bash
#
mkdir temp
cd temp
rm *
../f90split ../chrpak.f90
#
for FILE in `ls -1 *.f90`;
do
gfortran -c -g $FILE >& compiler.txt
if [ $? -ne 0 ]; then
echo "Errors compiling " $FILE
exit
fi
rm compiler.txt
done
rm *.f90
#
ar qc libchrpak.a *.o
rm *.o
#
mv libchrpak.a ../
#~/lib/$ARCH
cd ..
rmdir temp
#
echo "Library installed as ~/lib/$ARCH/libchrpak.a"
Libraries/f90split.f90 0 → 100644
View file @ cde68ad2
This diff is collapsed.
Libraries/f90split.sh 0 → 100644
View file @ cde68ad2
#!/bin/bash
#
gfortran -c -g f90split.f90 >& compiler.txt
if [ $? -ne 0 ]; then
echo "Errors compiling f90split.f90"
exit
fi
rm compiler.txt
#
gfortran f90split.o
if [ $? -ne 0 ]; then
echo "Errors linking and loading f90split.o"
exit
fi
rm f90split.o
#
chmod ugo+x a.out
mv a.out ~/bin/$ARCH/f90split
#
echo "Program installed as ~/bin/$ARCH/f90split"
Libraries/geometry.f90 0 → 100644
View file @ cde68ad2
This diff is collapsed.
Libraries/geometry.sh 0 → 100755
View file @ cde68ad2
#!/bin/bash
#
mkdir temp
cd temp
rm *
../f90split ../geometry.f90
#
for FILE in `ls -1 *.f90`;
do
gfortran -c -g $FILE >& compiler.txt
if [ $? -ne 0 ]; then
echo "Errors compiling " $FILE
exit
fi
rm compiler.txt
done
rm *.f90
#
ar qc libgeometry.a *.o
rm *.o
#
mv libgeometry.a ../
cd ..
rmdir temp
#
echo "Library installed as ~/lib/$ARCH/libgeometry.a"
README 0 → 100644
View file @ cde68ad2
#-------------------------------------------------------------------------------
# |
# POLYMER BUILDER INSTALLATION AND USAGE MANUAL |
# |
#-------------------------------------------------------------------------------
COMPILING
Two makefiles provided for gfortran and intel fortran compilers. These are
Makefile.gfortran : Makefile for gcc/gfortran compilers.
Makefile.ifort : Makefile for intel fortran compilers.
Enter the source directory:
cd Src
Copy appropriate make file to file name makefile:
cp Makefile.gfortran makefile
Clean the source directory:
make clean
Compile the source code:
make
This will create an executable file " polybuild " and it will copy this file to the test folder
cp polybuild ../Test
and copy polybuild to local bin folder if it exist:
cp polybuild ~/bin/
That completes building the source code.
USAGE:
Polybuild will need a configuration file which should be supplied via command line to execute:
polybuild config.cb
Configuration file contains setup information for the polybuild.
There are three modes of the program.
1- Polymer Building Mode
2- Composite Building Mode
3- Fiber Building Mode
These three modes can be combined except composite and fiber modes. Fiber modes can only be combined with Polymer mode.
1- Polymer Building Mode:
README.md 0 → 100644
View file @ cde68ad2
# PolyBuild
Polymer ceramic composites building tool for molecular simulations
Dundar Yilmaz
Penn State University
department of Mechanical and Nuclear Engineering
This tool creates complex structures for molecular simulations.
These structures are:
1- Amorphous polymers
2- Crystal structures
3- Fibers
4- Core-shell fibers
SRC_OLD/ceramic_type.f90 0 → 100644
View file @ cde68ad2
module ceramic_type
use checkbox
use IFPORT
implicit none
type ceramic
integer :: n_atoms
integer :: block_id
integer , allocatable,dimension(:) :: atype
integer , allocatable,dimension(:) :: deleted
integer , allocatable,dimension(:) :: surf_type
integer , allocatable,dimension(:) :: functional_site
character(len=3),allocatable , dimension (:) :: symbol
integer, allocatable, dimension(:) :: indx,chain
real*8 , allocatable , dimension(:,:) :: coord
real*8 , allocatable ,dimension(:) :: mass
real*8 , dimension(3,3) :: h_matrix
contains
procedure, public :: print_block
procedure , public :: create_surface_roughness
procedure , public :: find_functional_sites
procedure, public :: create_defects
end type
type(ceramic ) , allocatable , dimension(:) :: blocks
! type(ceramic ) :: base_ceramic
contains
function read_unit_cell(filename) result(new_ceramic)
type(ceramic) :: new_ceramic
character(len=*) :: filename
character(len=256) :: cwd,fullname
integer :: n_atoms,i
real*8 ,dimension(3):: vec
real*8 :: dumy
i=getcwd(cwd)
fullname=cwd(:LNBLNK(cwd))//"/"//filename(:LNBLNK(filename))
print*,fullname
open(13,file=filename)
read(13,*) n_atoms
new_ceramic%n_atoms=n_atoms
allocate(&
& new_ceramic%symbol(n_atoms) , &
& new_ceramic%indx(n_atoms) , &
& new_ceramic%mass(n_atoms) , &
& new_ceramic%coord(n_atoms,3), &
& new_ceramic%atype(n_atoms) ,&
& new_ceramic%chain(n_atoms) &
)
read(13,*) !scale
read(13,*) new_ceramic%h_matrix(1,1:3)
read(13,*) new_ceramic%h_matrix(2,1:3)
read(13,*) new_ceramic%h_matrix(3,1:3)
read(13,*) ! Direct
!print*,new_ceramic%h_matrix(1,1)
!print*,new_ceramic%h_matrix(2,2)
!print*,new_ceramic%h_matrix(3,3)
do i=1,n_atoms
read(13,*) new_ceramic%atype(i),vec,&
& new_ceramic%chain(i) , new_ceramic%symbol(i),new_ceramic%mass(i)
new_ceramic%coord(i,1) =vec(1)*new_ceramic%h_matrix(1,1)
new_ceramic%coord(i,2) =vec(2)*new_ceramic%h_matrix(2,2)
new_ceramic%coord(i,3) =vec(3)*new_ceramic%h_matrix(3,3)
end do
close(13)
end function
!=========================================================
function fiber_shell(center,radius_in,radius_out,base_ceramic) result( new_shell)
type(ceramic) :: new_shell,base_ceramic
real*8 ,dimension(2) :: center
real*8 , dimension(2,2):: bases
real*8 :: radius_out,dist,radius_in
integer :: nx,ny,nz
logical :: count_flag= .true.
integer :: cnt=0,ix,iy,iz,q,qb
real*8 , dimension(3) :: origin,vec
nz=int(BOX(3)/base_ceramic%h_matrix(3,3))
nx=int(radius_out/base_ceramic%h_matrix(1,1))
ny=int(radius_out/base_ceramic%h_matrix(2,2))
bases(1,1:2)=(/0.d0,0.d0/)
bases(2,1:2)=(/0.5d0,0.5d0/)
10 continue !print*,'fiber core creating '
do iz=1,nz
do ix=-nx,nx
do iy=-ny,ny
origin=(/center(1),center(2),0.d0/)+&
& (/ &
& dble(ix)*base_ceramic%h_matrix(1,1),&
& dble(iy)*base_ceramic%h_matrix(2,2),&
& dble(iz-1)*base_ceramic%h_matrix(3,3)/)
! base1
do qb=1,2
origin(1:2)=origin(1:2)+(/bases(qb,1)*base_ceramic%h_matrix(1,1),&
& bases(qb,2)*base_ceramic%h_matrix(2,2)/)
dist=sqrt(dot_product(origin(1:2)-center,origin(1:2)-center))
if ( dist > radius_out .or. dist < radius_in) cycle
!write(*,'(A,4f8.2)') 'Base : ',origin(1:2),dist
do q=1,base_ceramic%n_atoms
cnt=cnt+1
if( count_flag ) cycle
new_shell%symbol(cnt)=base_ceramic%symbol(q)
new_shell%atype(cnt) =base_ceramic%atype(q)
new_shell%mass(cnt) =base_ceramic%mass(q)
new_shell%coord(cnt,:)=base_ceramic%coord(q,:)+origin
end do
end do
end do
end do
end do
if( .not. count_flag) then
write(*,'(A,x,f6.2,x,A,x,f6.2)') 'Fiber shell created with inner radius ',radius_in, &
& 'outer radius',radius_out
return
end if
!print*,' Counted'
write(*,'(A,x,i8)')'Number of atoms at the shell region:', cnt
new_shell%n_atoms=cnt
allocate(new_shell%symbol(cnt) ,&
& new_shell%atype(cnt), &
& new_shell%mass(cnt), &
& new_shell%deleted(cnt), &
& new_shell%indx(cnt), &
& new_shell%functional_site(cnt), &
& new_shell%coord(cnt,3))
cnt=0
count_flag=.false.
goto 10
end function
!===================================================
function fiber_core(center,radius,base_ceramic) result( new_core)
type(ceramic) :: new_core,base_ceramic
real*8 ,dimension(2) :: center
real*8 , dimension(2,2):: bases
real*8 :: radius,dist
integer :: nx,ny,nz
logical :: count_flag= .true.
integer :: cnt=0,ix,iy,iz,q,qb
real*8 , dimension(3) :: origin,vec
nz=int(BOX(3)/base_ceramic%h_matrix(3,3))
nx=int(radius/base_ceramic%h_matrix(1,1))
ny=int(radius/base_ceramic%h_matrix(2,2))
bases(1,1:2)=(/0.d0,0.d0/)
bases(2,1:2)=(/0.5d0,0.5d0/)
10 continue !print*,'fiber core creating '
do iz=1,nz
do ix=-nx,nx
do iy=-ny,ny
origin=(/center(1),center(2),0.d0/)+&
& (/ &
& dble(ix)*base_ceramic%h_matrix(1,1),&
& dble(iy)*base_ceramic%h_matrix(2,2),&
& dble(iz-1)*base_ceramic%h_matrix(3,3)/)
! base1
do qb=1,2
origin(1:2)=origin(1:2)+(/bases(qb,1)*base_ceramic%h_matrix(1,1),&
& bases(qb,2)*base_ceramic%h_matrix(2,2)/)
dist=sqrt(dot_product(origin(1:2)-center,origin(1:2)-center))
if ( dist > radius) cycle
!write(*,'(A,4f8.2)') 'Base : ',origin(1:2),dist
do q=1,base_ceramic%n_atoms
cnt=cnt+1
if( count_flag ) cycle
new_core%symbol(cnt)=base_ceramic%symbol(q)
new_core%atype(cnt) =base_ceramic%atype(q)
new_core%mass(cnt) =base_ceramic%mass(q)
new_core%coord(cnt,:)=base_ceramic%coord(q,:)+origin
end do
end do
end do
end do
end do
if( .not. count_flag) then
write(*,'(A)') 'Fiber core created'
return
end if
!print*,' Counted'
write(*,'(A,i8)')'Total number of atoms : ', cnt
new_core%n_atoms=cnt
allocate(new_core%symbol(cnt) ,&
& new_core%atype(cnt), &
& new_core%mass(cnt), &
& new_core%deleted(cnt), &
& new_core%indx(cnt), &
& new_core%functional_site(cnt), &
& new_core%coord(cnt,3))
cnt=0
count_flag=.false.
goto 10
end function
function fill_blocks(box_id,base_ceramic) result(new_block)
type(ceramic) :: new_block,base_ceramic
integer :: box_id
integer :: nx,ny,nz
integer :: i , j ,k ,q
integer :: n_atoms,n_ucell,cnt
real*8 :: r_nx, r_ny, r_nz
real*8 , dimension(3) :: pos
real*8 , dimension(3) :: vec
real*8 , dimension(3,3) :: h_matrix
real*8 :: surf_depth=2.d0
real*8 :: z_cor, z_limit
nx=nint(box_dimensions(box_id,4)/base_ceramic%h_matrix(1,1))
ny=nint(box_dimensions(box_id,5)/base_ceramic%h_matrix(2,2))
nz=nint(box_dimensions(box_id,6)/base_ceramic%h_matrix(3,3))
r_nx=box_dimensions(box_id,4)/base_ceramic%h_matrix(1,1)
r_ny=box_dimensions(box_id,5)/base_ceramic%h_matrix(2,2)
r_nz=box_dimensions(box_id,6)/base_ceramic%h_matrix(3,3)
h_matrix=base_ceramic%h_matrix
print*, box_id
print*,'X: ', box_dimensions(box_id,4)/base_ceramic%h_matrix(1,1)
print*,'Y :', box_dimensions(box_id,5)/base_ceramic%h_matrix(2,2)
print*, 'Z' , box_dimensions(box_id,6)/base_ceramic%h_matrix(3,3)
write(*,*) nx,ny,nz
write(*,*) nx*ny*nz*base_ceramic%n_atoms
if ( abs(r_nx-nx) .gt. 0.01) then
write(*,'(A,x,f8.4,x,i4)') 'Lattice mismatch in X direction : ', r_nx,nx
h_matrix(1,1)=box_dimensions(box_id,4)/dble(nx)
write(*,'(A,x,f8.4,x,A,x,f8.4,x,A,x,f5.2)') &
& 'Old lattice constant in x:', base_ceramic%h_matrix(1,1),&
& 'New Lattice Constant :',h_matrix(1,1),&
& 'Strain applied:(%)',(base_ceramic%h_matrix(1,1)-h_matrix(1,1))/h_matrix(1,1)*100
end if
if ( abs(r_ny-ny) .gt. 0.01) then
write(*,'(A,x,f8.4,x,i4)') 'Lattice mismatch in Y direction : ', r_ny,ny
h_matrix(2,2)=box_dimensions(box_id,5)/dble(ny)
write(*,'(A,x,f8.4,x,A,x,f8.4,x,A,x,f5.2)') &
& 'Old lattice constant in y:', base_ceramic%h_matrix(2,2),&
& 'New Lattice Constant :',h_matrix(2,2),&
& 'Strain applied:(%)',(base_ceramic%h_matrix(2,2)-h_matrix(2,2))/h_matrix(2,2)*100
end if
if ( abs(r_nz-nz) .gt. 0.01) then
write(*,'(A,x,f8.4,x,i4)') 'Lattice mismatch in Z direction : ', r_nz,nz
h_matrix(3,3)=box_dimensions(box_id,6)/dble(nz)
write(*,'(A,x,f8.4,x,A,x,f8.4,x,A,x,f5.2)') &
& 'Old lattice constant in z:', base_ceramic%h_matrix(3,3),&
& 'New Lattice Constant :',h_matrix(3,3),&
& 'Strain applied:(%)',(base_ceramic%h_matrix(3,3)-h_matrix(3,3))/h_matrix(3,3)*100
end if
! nx=int(box_lengths(1)/base_ceramic%h_matrix(1,1))
! ny=int(box_lengths(2)/base_ceramic%h_matrix(2,2))
! nz=int(box_lengths(3)/base_ceramic%h_matrix(3,3))
n_atoms=nx*ny*nz*base_ceramic%n_atoms
n_ucell=base_ceramic%n_atoms
new_block%n_atoms=n_atoms
new_block%block_id=box_id
write(*,'(A,x,i2)')'Filling Block : ', box_id
!print*,' nx : ',nx,' ny : ',ny,' nz : ',nz
!print*,(box_lengths(1)/base_ceramic%h_matrix(1,1))
!print*,' Number of atoms in block : ',n_atoms
!print*,'Size-x : ',dble(nx)*base_ceramic%h_matrix(1,1)
!print*,'Size-y : ',dble(ny)*base_ceramic%h_matrix(2,2)
!print*,'Size-z : ',dble(nz)*base_ceramic%h_matrix(3,3)
allocate(&
& new_block%symbol(n_atoms) , &
& new_block%indx(n_atoms) , &
& new_block%mass(n_atoms) , &
& new_block%coord(n_atoms,3), &
& new_block%atype(n_atoms), &
& new_block%deleted(n_atoms), &
& new_block%surf_type(n_atoms), &
& new_block%functional_site(n_atoms) &
)
new_block%surf_type(:)=0
new_block%functional_site(:)=0
new_block%deleted(:)=0
cnt=0
do i=1,nx
do j=1,ny
do k=1,nz
! print*,i,j,k
pos=box_dimensions(box_id,:)+&
(/&
& dble(i-1)*h_matrix(1,1),&
& dble(j-1)*h_matrix(2,2),&
& dble(k-1)*h_matrix(3,3)/)
! print*,i , j ,k,box_id
! print*,box_dimensions(box_id,:)
! print*,pos
do q=1,base_ceramic%n_atoms
! if( (i .eq. 1) .and. (base_ceramic%chain(q) .eq. 1) ) cycle
cnt=cnt+1
new_block%symbol(cnt)=base_ceramic%symbol(q)
new_block%atype(cnt) =base_ceramic%atype(q)
new_block%mass(cnt) =base_ceramic%mass(q)
new_block%coord(cnt,:)=base_ceramic%coord(q,:)+pos
! new_block%coord((cnt-1)*n_ucell+1:cnt*n_ucell,1)=base_ceramic%coord(:,1)+pos(1)
! new_block%coord((cnt-1)*n_ucell+1:cnt*n_ucell,2)=base_ceramic%coord(:,2)+pos(2)
! new_block%coord((cnt-1)*n_ucell+1:cnt*n_ucell,3)=base_ceramic%coord(:,3)+pos(3)
end do
end do
end do
end do
n_atoms=cnt
!print*, 'N atoms : ',cnt
! pause
!detect surface atoms
!-----------------------------------------------
! surf_type = 1 upper xy plane
! surf_type = 2 lower xy plane
! surf_type = 0 bulk
!-----------------------------------------------
do i=1,n_atoms
vec= new_block%coord(i,:)
z_cor=new_block%coord(i,3)
z_limit=box_dimensions(new_block%block_id,3)+box_lengths(3)-surf_depth
if(z_cor > z_limit) new_block%surf_type(i)=1
z_limit=box_dimensions(new_block%block_id,3)+surf_depth
if(z_cor < z_limit) new_block%surf_type(i)=2
end do
end function
function print_block(this,filename,indx) result(iostat)
class(ceramic) :: this
character(len=*) :: filename
integer :: indx,cnt,i
integer :: iostat
if( indx == 0) then
!---- print header ---- '