add angle force in polymer

examples/USER/sph/sdpd-polymer-3D-test/run3d.sh

@@ -7,7 +7,7 @@
 ../../../../tools/restart2data poly3d.restart poly3d.txt
 
 
- awk -v cutoff=3.0 -v Nbeads=1 -v Nsolvent=1 -v Npoly=full \
+ awk -v cutoff=3.0 -v Nbeads=6 -v Nsolvent=2 -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

examples/USER/sph/sdpd-polymer-3D-test/sdpd-polymer3D-inti.lmp

@@ -13,8 +13,8 @@ region             box block 0.0 ${L}  0.0 ${L} 0.0 ${L} units box
 create_box         1 box
 
 # create fluid particles
-lattice            fcc 0.025e-3
-#lattice            sc 0.025e-3
+#lattice            fcc 0.025e-3
+lattice            sc 0.025e-3
 create_atoms       1 region box
 
 variable           sdpd_mass equal 3.90625e-12

examples/USER/sph/sdpd-polymer-3D-test/sdpd-polymer3D-run.lmp

@@ -15,9 +15,9 @@ read_data           poly3d.txt
 lattice            sc 0.025e-3
 create_atoms       1 region box
 
-#special_bonds   lj 1 1 1
-#bond_style fene
-#bond_coeff * 3e-2 5e-5 0 0
+special_bonds   lj 1 1 1
+bond_style fene
+bond_coeff * 1 5e-5 0 0
 
 
 mass               1 3.90625e-12
@@ -31,9 +31,10 @@ variable          h equal 6.5e-5
 # 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
 
-variable           sdpd_temp index 1e1
-pair_style         hybrid sdpd
-pair_coeff         * *    sdpd 1000 0.1 5.0e-3 6.5e-5 ${sdpd_temp}
+variable           sdpd_temp index 1e-9
+pair_style         hybrid  sdpd
+#pair_coeff         * * sph/rhosum   ${h}
+pair_coeff         * *    sdpd 1000 1e-4 1.0e-4 6.5e-5 ${sdpd_temp}
 
 compute            rho_peratom all meso_rho/atom
 compute            e_peratom all meso_e/atom
@@ -59,5 +60,5 @@ thermo             1
 
 
 neighbor           3.0e-6 bin
-timestep           1e-5
+timestep           1e-7
 run                40000

examples/USER/sph/sdpd-polymer-angle/addpolymer.awk

@@ -0,0 +1,128 @@
+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"
+  print "1 angle types"
+  next
+}
+
+/atoms/{
+  natoms=$1
+  print
+  printf("%s bonds\n", "_NUMBER_OF_BOUNDS_")
+  printf("%s angles\n", "_NUMBER_OF_ANGLES_")
+
+  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
+    } 
+ }
+
+ printf("\nAngles\n\n")
+
+ iangle = 0
+
+  for (a=1; a<iatom; a++) {
+    if (isbound(a+1)&&isbound(a)) {
+     iangle++
+      ia = a
+      ja = a+1
+      ka = a+2
+      angletype=1
+      print iangle, angletype, ia, ja, ka
+    }
+  }
+
+}