Permalink
Browse files

Create a template for wall bounded polymer solution

  • Loading branch information...
1 parent 211537d commit 6fc819f99a702f3ae16c3f2d6cf73d483df27595 Xie Qingguang committed Dec 16, 2011
View
109 examples/USER/sph/sdpd-polymer-wall/addpolymer.awk
@@ -0,0 +1,109 @@
+function fabs(var) {
+ return var>0?var:-var
+}
+
+function isbound(atom_number, period, rem, current_npoly) {
+ period = Nbeads + Nsolvent
+ rem = (atom_number-1)%(period) # from 0 to period-1
+ current_npoly = int(atom_number/period) + 1
+ return (rem<Nbeads-1) && (atom_number<iatom) && (current_npoly<=Npoly)
+}
+
+BEGIN {
+ inatoms=0
+ lo=1; hi=2
+ x=1; y=2; z=3
+ iatom=0
+ if (Npoly=="full") {
+ Npoly = 1e22
+ }
+}
+
+/LAMMPS/{
+ print
+ next
+}
+
+/xlo xhi/{
+ box[x,lo]=$1
+ box[x,hi]=$1
+}
+
+/ylo yhi/{
+ box[y,lo]=$1
+ box[y,hi]=$2
+}
+
+/zlo zhi/{
+ box[z,lo]=$1
+ box[z,hi]=$2
+}
+
+/atom types/{
+ print
+ print "1 bond types"
+ next
+}
+
+/atoms/{
+ natoms=$1
+ print
+ printf("%s bonds\n", "_NUMBER_OF_BOUNDS_")
+ next
+}
+
+(NF>0)&&($1=="Atoms"){
+ inatoms=1
+ print
+ # skip empty line
+ getline
+ printf "\n"
+ next
+}
+
+inatoms && (NF==0) {
+ inatoms = 0
+ print
+ next
+}
+
+inatoms{
+ # here I get one atom
+ iatom++
+ R[x]=$4; R[y]=$5; R[z]=$6
+ if (iatom>1) {
+ for (idim=1; idim<=3; idim++) {
+ if (fabs(R[idim]- prevR[idim])>cutoff) {
+ if (R[idim]<prevR[idim]) image[idim]++; else image[idim]--
+ }
+ }
+ } else {
+ image[x]=0; image[y]=0; image[z]=0
+ }
+ prevR[x]=R[x]; prevR[y]=R[y]; prevR[z]=R[z]
+ # change image field
+ $(NF-2)=image[x]; $(NF-1)=image[y]; $(NF)=image[z];
+ # add molecule ID
+ # $6=$6 " 0"
+ print $0
+ next
+}
+
+!inatoms {
+ print
+ next
+}
+
+END {
+ printf("\nBonds\n\n")
+ ibond = 0
+ for (q=1; q<iatom; q++) {
+ if (isbound(q)) {
+ ibond++
+ ip = q
+ jp = q+1
+ bondtype=1
+ print ibond, bondtype, ip, jp
+ }
+ }
+}
View
3 examples/USER/sph/sdpd-polymer-wall/plotpunto.sh
@@ -0,0 +1,3 @@
+bash post.sh
+punto -s 2 -G [-1e-3:1e-3] -D 2 -c 3 punto.dat
+
View
6 examples/USER/sph/sdpd-polymer-wall/post.sh
@@ -0,0 +1,6 @@
+#! /bin/bash
+rm -rf punto.dat
+#awk 'fl{print $3, $4, $5, $6} /ITEM: ATOMS/{fl=1}' *.dat > punto.dat
+for filename in dump*.dat; do echo $filename; sed '1,8d' $filename > ${filename/.dat/.del};done
+awk 'fl{print $1, $2, $3, $4} /ITEM: ATOMS/{fl=1}' *.del > punto.dat
+
View
17 examples/USER/sph/sdpd-polymer-wall/run.sh
@@ -0,0 +1,17 @@
+#! /bin/bash
+
+
+rm -rf dump*
+rm -rf image*
+
+../../../../src/lmp_linux -in sdpd-polymer-inti.lmp
+../../../../tools/restart2data poly.restart poly.txt
+
+
+ awk -v cutoff=3.0 -v Nbeads=18 -v Nsolvent=18 -v Npoly=full \
+ -f addpolymer.awk poly.txt > poly2.txt
+ nbound=$(tail -n 1 poly2.txt | awk '{print $1}')
+ sed "s/_NUMBER_OF_BOUNDS_/$nbound/1" poly2.txt > poly.txt
+
+time /scratch/qingguang/prefix-nana/bin/mpirun -np 6 ../../../../src/lmp_linux -in sdpd-polymer-run.lmp
+#../../../../src/lmp_linux -in sdpd-polymer-run.lmp
View
16 examples/USER/sph/sdpd-polymer-wall/run3d.sh
@@ -0,0 +1,16 @@
+#! /bin/bash
+
+
+
+
+../../../../src/lmp_linux -in sdpd-polymer3D-inti.lmp
+../../../../tools/restart2data poly3d.restart poly3d.txt
+
+
+ awk -v cutoff=3.0 -v Nbeads=3 -v Nsolvent=1 -v Npoly=full \
+ -f addpolymer.awk poly3d.txt > poly3.txt
+ nbound=$(tail -n 1 poly3.txt | awk '{print $1}')
+ sed "s/_NUMBER_OF_BOUNDS_/$nbound/1" poly3.txt > poly3d.txt
+
+mpirun -np 2 ../../../../src/lmp_linux -in sdpd-polymer3D-run.lmp
+#../../../../src/lmp_linux -in sdpd-polymer-run.lmp
View
65 examples/USER/sph/sdpd-polymer-wall/sdpd-polymer-inti.lmp
@@ -0,0 +1,65 @@
+echo both
+#atom_modify sort 3.0e-6 3.0e-6
+dimension 2
+units si
+atom_style hybrid bond meso
+boundary p p p
+
+# create simulation box
+#2D box
+variable Lx equal 2.000e-3
+region box block 0.0000e-3 ${Lx} 0.0000e-3 ${Lx} -1.0e-6 1.0e-6 units box
+create_box 1 box
+
+# create fluid particles
+lattice sq 0.025e-3
+create_atoms 1 region box
+
+#set bead number of polymers and free particles
+#attention the relation between sum(Sum) of beadnum with free and Lx/latticelength(N) satisfy N%Sum=0
+#beadnum 5
+#freenum 5
+#string coefficient H
+#Maxiaml distance R0
+
+variable sdpd_mass equal 2.776E-6
+mass 1 ${sdpd_mass}
+#set group all meso_rho 1000.0
+
+# use Tait's EOS in combination with Morris' laminar viscosity.
+# We set rho_0 = 1000 kg/m^3, c = 0.1 m/s, h = 6.5e-5 m.
+# The dynamic viscosity is set to 1.0e-3 Pa s, corresponding to a kinematic viscosity of 1.0e-6 m^2/s
+#The temperature is set to 1e5C.
+variable sdpd_temp index 0.0
+#pair_style hybrid sdpd
+#pair_coeff * * sdpd 1000 0.1 5.0e-3 6.5e-5 ${sdpd_temp}
+
+# compute rho_peratom all meso_rho/atom
+# compute e_peratom all meso_e/atom
+# compute ke_peratom all ke/atom
+# compute esph all reduce sum c_e_peratom
+# compute ke all ke
+# compute sdpd_kin all temp
+# variable etot equal c_ke+c_esph
+
+
+# do full time integration for shear driver and fluid, but keep walls stationary
+#fix integrate_fix_full all meso
+#variable A_kol equal 1.0
+#variable var_kol atom ${A_kol}*${sdpd_mass}*sin(2*PI*y/${Lx})
+#fix kol_force all addforce v_var_kol 0.0 0.0
+
+
+# dump dump_id all custom 100 dump*.dat id type x y vx vy c_rho_peratom c_e_peratom
+# dump_modify dump_id first yes
+# dump_modify dump_id sort id
+# dump_modify dump_id pad 8
+# thermo 10
+# thermo_style custom step c_sdpd_kin
+# thermo_modify norm no
+
+#neighbor 1.0e-6 bin
+timestep 0.0e-5
+write_restart poly.restart
+
+#run 0
View
71 examples/USER/sph/sdpd-polymer-wall/sdpd-polymer-run.lmp
@@ -0,0 +1,71 @@
+echo both
+#atom_modify sort 3.0e-6 3.0e-6
+dimension 2
+units si
+atom_style hybrid bond meso
+boundary p p p
+
+# create simulation box
+#2D box
+variable Lx equal 2.00e-3
+
+region box block 0.0000e-3 ${Lx} 0.0000e-3 ${Lx} -1.0e-6 1.0e-6 units box
+read_data poly.txt
+
+# create fluid particles
+lattice sq 0.025e-3
+bond_style fene
+special_bonds lj 1 1 1i
+#H*ro^2/(kT)=100
+variable H_fene equal 5
+#variable H_fene equal 0.5e-1
+variable r0_fene equal 5e-5
+bond_coeff 1 ${H_fene} ${r0_fene} 0.0 0.0
+
+
+variable sdpd_mass equal 2.778E-6
+mass 1 ${sdpd_mass}
+set group all meso_rho 1000.0
+
+# use Tait's EOS in combination with Morris' laminar viscosity.
+# We set rho_0 = 1000 kg/m^3, c = 0.1 m/s, h = 6.5e-5 m.
+# The dynamic viscosity is set to 1.0e-3 Pa s, corresponding to a kinematic viscosity of 1.0e-6 m^2/s
+#The temperature is set to 1e14.
+variable sdpd_temp index 1e13
+pair_style hybrid sdpd
+pair_coeff * * sdpd 1000 0.1 1.0e-3 6.5e-5 ${sdpd_temp}
+
+compute rho_peratom all meso_rho/atom
+compute e_peratom all meso_e/atom
+compute ke_peratom all ke/atom
+compute esph all reduce sum c_e_peratom
+compute ke all ke
+compute sdpd_kin all temp
+variable etot equal c_ke+c_esph
+
+
+# do full time integration for shear driver and fluid, but keep walls stationary
+fix integrate_fix_full all meso
+variable A_kol equal 0.1
+variable var_kol atom ${A_kol}*${sdpd_mass}*cos(2*PI*y/${Lx})
+fix kol_force all addforce v_var_kol 0.0 0.0
+
+
+#dump dump_id all custom 10000 dump*.dat id type x y vx vy c_rho_peratom c_e_peratom
+dump dump_id all custom 1000 dump*.dat x y vx vy id type c_rho_peratom
+dump_modify dump_id first yes
+dump_modify dump_id sort id
+dump_modify dump_id pad 8
+thermo 10
+thermo_style custom step c_sdpd_kin
+thermo_modify norm no
+
+dump imgDump all image 1000 image.*.jpg type type bond atom 1e-5 adiam 1e-5 view 0 0 zoom 1.75
+dump_modify imgDump pad 9
+
+neighbor 1.0e-6 bin
+#timestep <0.5*delta/Cs
+#include detime
+timestep 1.0e-5
+
+run 1000000000
View
63 examples/USER/sph/sdpd-polymer-wall/sdpd-polymer.lmp
@@ -0,0 +1,63 @@
+echo both
+dimension 2
+units si
+atom_style meso
+boundary p p p
+
+# create simulation box
+#2D box
+variable Lx equal 1.000e-3
+region box block 0.0000e-3 ${Lx} 0.0000e-3 ${Lx} -1.0e-6 1.0e-6 units box
+create_box 1 box
+
+# create fluid particles
+lattice sq 0.025e-3
+create_atoms 1 region box
+
+#set bead number of polymers and free particles
+#attention the relation between sum(Sum) of beadnum with free and Lx/latticelength(N) satisfy N%Sum=0
+#beadnum 5
+#freenum 5
+#string coefficient H
+#Maxiaml distance R0
+
+variable sdpd_mass equal 6.25E-7
+mass 1 ${sdpd_mass}
+set group all meso_rho 1000.0
+
+# use Tait's EOS in combination with Morris' laminar viscosity.
+# We set rho_0 = 1000 kg/m^3, c = 0.1 m/s, h = 6.5e-5 m.
+# The dynamic viscosity is set to 1.0e-3 Pa s, corresponding to a kinematic viscosity of 1.0e-6 m^2/s
+#The temperature is set to 1e5C.
+variable sdpd_temp index 0.0
+pair_style hybrid sdpd
+pair_coeff * * sdpd 1000 0.1 5.0e-3 6.5e-5 ${sdpd_temp}
+
+compute rho_peratom all meso_rho/atom
+compute e_peratom all meso_e/atom
+compute ke_peratom all ke/atom
+compute esph all reduce sum c_e_peratom
+compute ke all ke
+compute sdpd_kin all temp
+variable etot equal c_ke+c_esph
+
+
+# do full time integration for shear driver and fluid, but keep walls stationary
+fix integrate_fix_full all meso
+#variable A_kol equal 1.0
+#variable var_kol atom ${A_kol}*${sdpd_mass}*sin(2*PI*y/${Lx})
+#fix kol_force all addforce v_var_kol 0.0 0.0
+
+
+dump dump_id all custom 100 dump*.dat id type x y vx vy c_rho_peratom c_e_peratom
+dump_modify dump_id first yes
+dump_modify dump_id sort id
+dump_modify dump_id pad 8
+thermo 10
+thermo_style custom step c_sdpd_kin
+thermo_modify norm no
+
+neighbor 3.0e-6 bin
+timestep 8.0e-5
+
+run 40000 every 1000 "velocity all zero linear"
View
7 examples/USER/sph/sdpd-polymer/run.sh
@@ -1,16 +1,17 @@
#! /bin/bash
-
+rm -rf dump*
+rm -rf image*
../../../../src/lmp_linux -in sdpd-polymer-inti.lmp
../../../../tools/restart2data poly.restart poly.txt
- awk -v cutoff=3.0 -v Nbeads=9 -v Nsolvent=18 -v Npoly=full \
+ awk -v cutoff=3.0 -v Nbeads=18 -v Nsolvent=18 -v Npoly=full \
-f addpolymer.awk poly.txt > poly2.txt
nbound=$(tail -n 1 poly2.txt | awk '{print $1}')
sed "s/_NUMBER_OF_BOUNDS_/$nbound/1" poly2.txt > poly.txt
-time /scratch/qingguang/prefix-nana/bin/mpirun -np 4 ../../../../src/lmp_linux -in sdpd-polymer-run.lmp
+time /scratch/qingguang/prefix-nana/bin/mpirun -np 6 ../../../../src/lmp_linux -in sdpd-polymer-run.lmp
#../../../../src/lmp_linux -in sdpd-polymer-run.lmp
View
2 examples/USER/sph/sdpd-polymer/sdpd-polymer-inti.lmp
@@ -7,7 +7,7 @@ boundary p p p
# create simulation box
#2D box
-variable Lx equal 1.000e-3
+variable Lx equal 2.000e-3
region box block 0.0000e-3 ${Lx} 0.0000e-3 ${Lx} -1.0e-6 1.0e-6 units box
create_box 1 box
View
26 examples/USER/sph/sdpd-polymer/sdpd-polymer-run.lmp
@@ -7,16 +7,18 @@ boundary p p p
# create simulation box
#2D box
-variable Lx equal 1.00e-3
+variable Lx equal 2.00e-3
region box block 0.0000e-3 ${Lx} 0.0000e-3 ${Lx} -1.0e-6 1.0e-6 units box
read_data poly.txt
# create fluid particles
lattice sq 0.025e-3
bond_style fene
-special_bonds lj 1 1 1
-variable H_fene equal 1e-1
+special_bonds lj 1 1 1i
+#H*ro^2/(kT)=100
+variable H_fene equal 5
+#variable H_fene equal 0.5e-1
variable r0_fene equal 5e-5
bond_coeff 1 ${H_fene} ${r0_fene} 0.0 0.0
@@ -28,8 +30,8 @@ set group all meso_rho 1000.0
# use Tait's EOS in combination with Morris' laminar viscosity.
# We set rho_0 = 1000 kg/m^3, c = 0.1 m/s, h = 6.5e-5 m.
# The dynamic viscosity is set to 1.0e-3 Pa s, corresponding to a kinematic viscosity of 1.0e-6 m^2/s
-#The temperature is set to 1e5C.
-variable sdpd_temp index 1e12
+#The temperature is set to 1e14.
+variable sdpd_temp index 1e13
pair_style hybrid sdpd
pair_coeff * * sdpd 1000 0.1 1.0e-3 6.5e-5 ${sdpd_temp}
@@ -44,24 +46,26 @@ variable etot equal c_ke+c_esph
# do full time integration for shear driver and fluid, but keep walls stationary
fix integrate_fix_full all meso
-variable A_kol equal 1
+variable A_kol equal 0.1
variable var_kol atom ${A_kol}*${sdpd_mass}*cos(2*PI*y/${Lx})
fix kol_force all addforce v_var_kol 0.0 0.0
#dump dump_id all custom 10000 dump*.dat id type x y vx vy c_rho_peratom c_e_peratom
-dump dump_id all custom 10000 dump*.dat x y vx vy id type c_rho_peratom
+dump dump_id all custom 1000 dump*.dat x y vx vy id type c_rho_peratom
dump_modify dump_id first yes
dump_modify dump_id sort id
dump_modify dump_id pad 8
thermo 10
thermo_style custom step c_sdpd_kin
thermo_modify norm no
-dump imgDump all image 10000 image.*.jpg type type bond atom 1e-5 adiam 1e-5 view 0 0 zoom 1.75
-dump_modify imgDump pad 8
+dump imgDump all image 1000 image.*.jpg type type bond atom 1e-5 adiam 1e-5 view 0 0 zoom 1.75
+dump_modify imgDump pad 9
neighbor 1.0e-6 bin
-timestep 1.0e-6
+#timestep <0.5*delta/Cs
+#include detime
+timestep 1.0e-5
-run 1000000
+run 1000000000

0 comments on commit 6fc819f

Please sign in to comment.