Permalink
Browse files

fix temperature

  • Loading branch information...
1 parent 286433a commit 7b529dcbb271205a1bfad230b9bbbd36280e0c8a Xie Qingguang committed Feb 14, 2012
@@ -1,56 +1,59 @@
+
function ext()
% get polymer configurations
- flist=dir("pdata/poly.1");
+ flist=dir("poly.1");
%nfile is the number of polymers
nfile = size(flist,1);
%Loop for all polymers
for kk=1:nfile
- % for kk=1:1
- fname=fullfile("pdata", flist(kk).name);
-%only one polymer
- if (kk==1)
- extfun = getpolymerext(fname);
-
-%multi polymers,then do sum
- else
- printf("kk=%d\n", kk);
- extfun = extfun + getpolymerext(fname);
+ fname=fullfile(".", flist(kk).name);
+ data = file2data(fname);
+ if (~exist("extfun", "var"))
+ extfun = getpolymerext(data);
+ rg2 = getrg2(data);
+ else
+ extfun = extfun + getpolymerext(data);
+ rg2 = rg2 + getpolymerext(data);
endif
- end
+ endfor
+
nb=getnbead(fname);
printf("Polymer with beads %d",nb);
printf(" the end-to-end distance is:");
-extfun=extfun/nfile
- % dtime = 0:size(extfun, 1)-1;
- %dlmwrite( "extx.dat", [dtime', extfun/nfile], ' ', "precision", "%e");
-
+extfun=extfun/nfile;
+rg2 = rg2/nfile;
+dtime = 0:size(extfun, 1)-1;
+ dlmwrite( "extx.dat", [dtime', extfun], ' ', "precision", "%e");
+ dlmwrite( "rg2.dat", [dtime', rg2], ' ', "precision", "%e");
endfunction
-function [extfun] = getpolymerext(fname)
+function data = file2data(fname)
+ if (~exist(fname, "file"))
+ error("cannot one file: %s\n", fname);
+ endif
% get number of beads
nb=getnbead(fname);
+ warning("ID1", "number of beads is %d", nb);
dim=3;
% read data without space between rows
data = dlmread(fname);
% keep only x and y and z
data = data(:, 1:dim);
-%put different snapshot in cells the order is column order 1xyz2xyz3xyz.........
-%the endend value is endbead-z
nsnap = size(data, 1)/nb;
data = permute(reshape(data(:,:)', 3, nb, nsnap), [3, 2, 1]);
+ warning("data structure: [%d %d %d]", size(data));
+ % data structure is
+ % data ( isnaphot, ibead, dim )
+endfunction
- % rdata structure is
- % rdata ( isnaphot, ibead, dim )
-
-
+function [extfun] = getpolymerext(data)
% the number of snapshots
%the end-to-end distance
-
Rhead = squeeze(data(:, 1, :));
Rtail = squeeze(data(:, end, :));
Re2e = Rtail - Rhead;
@@ -64,15 +67,21 @@ function ext()
Rg2 = mean(sumsq(Rincm, 3), 2);
endfunction
+function Rext = Rext(data, dim)
+ nb = size(data, dim);
+ R = squeeze(data(:, :, dim));
+
+ Rext = max(R') - min(R');
+endfunction
+
% get the number of beads
function [nbead] = getnbead(fname)
-fid = fopen (fname);
+ fid = fopen (fname);
nbead=-1;
do
line = strtrim(fgetl(fid));
nbead=nbead+1;
until strcmp(line, "")
fclose (fid);
-%printf("Polymer with ",nbead,"beads");
endfunction
@@ -1,5 +0,0 @@
-#! /bin/bash
-
-
-
-../../../../src/lmp_linux -in sdpd-polymer3D-run.lmp
@@ -1,40 +0,0 @@
-echo both
-dimension 3
-units si
-atom_style meso
-boundary f f f
-
-# create simulation box
-#3D box
-variable L equal 0.1e-4
-region box block 0.0 ${L} 0.0 ${L} 0.0 ${L} units box
-create_box 1 box
-
-mass 1 3.90625e-12
-set group all meso_rho 1000.0
-
-# smothing length
-variable h equal 6.5e-5
-
-
-# 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
-
-variable sdpd_temp index 1e-9
-pair_style sdpd
-pair_coeff * * 1000 1e-4 1.0e-4 6.5e-5 ${sdpd_temp}
-
-compute rho_peratom all meso_rho/atom
-compute sdpd_kin all temp
-
-# do full time integration for shear driver and fluid, but keep walls stationary
-fix integrate_fix_full all meso
-
-thermo_style custom step c_sdpd_kin
-thermo_modify norm no
-thermo 1
-
-neighbor 3.0e-6 bin
-timestep 0
-run 1
View
No changes.
View
No changes.
@@ -10,7 +10,7 @@ boundary p p p
# create simulation box
#3D box
variable L equal 3.0e-4
-region box block 0.0 ${L} 0.0 ${L} 0.0 ${L} units box
+region box prism 0.0 ${L} 0.0 ${L} 0.0 ${L} 0 0 0 units box
create_box 1 box
# create fluid particles
@@ -6,7 +6,6 @@ boundary p p p
# create simulation box
#3D box
variable L equal 3e-4
-region box block 0.0 ${L} 0.0 ${L} 0.0 ${L} units box
read_data poly3d.txt
special_bonds lj 1 1 1
@@ -37,12 +36,14 @@ compute ke all ke
compute sdpd_kin all temp
variable etot equal c_ke+c_esph
+# Lees Edwards type boundary conditions
+fix 2 all deform 1 xy erate 1000000 remap v
# do full time integration for shear driver and fluid, but keep walls stationary
fix integrate_fix_full all meso
-variable A_kol equal 5
-variable var_kol atom ${A_kol}*${sdpd_mass}*cos(2*PI*y/${L})
-fix kol_force all addforce v_var_kol 0.0 0.0
+# variable A_kol equal 5
+# variable var_kol atom ${A_kol}*${sdpd_mass}*cos(2*PI*y/${L})
+# fix kol_force all addforce v_var_kol 0.0 0.0
dump dump_id all custom 100000 dump*.dat id type x y z vx vy vz c_rho_peratom
@@ -52,7 +53,7 @@ dump_modify dump_id pad 8
thermo_style custom step c_sdpd_kin
thermo_modify norm no
thermo 1
-dump imgDump all image 100000 image.*.jpg type type atom no bond atom 5e-6 adiam 1e-5 view 60 60 zoom 1.3 box yes 0.01
+dump imgDump all image 100 image.*.jpg type type atom yes bond atom 5e-6 adiam 1e-5 view 60 60 zoom 1.3 box yes 0.01
dump_modify imgDump pad 9
@@ -1,13 +1,13 @@
function ext()
% get polymer configurations
- flist=dir("pdata/poly.*");
+ flist=dir("pdata/poly.1");
%nfile is the number of polymers
nfile = size(flist,1);
%Loop for all polymers
- for kk=1:nfile-1
+ for kk=1:nfile
% for kk=1:1
fname=fullfile("pdata", flist(kk).name);
%only one polymer
@@ -34,34 +34,37 @@ function ext()
nb=getnbead(fname);
dim=3;
% read data without space between rows
+
data = dlmread(fname);
% keep only x and y and z
data = data(:, 1:dim);
%put different snapshot in cells the order is column order 1xyz2xyz3xyz.........
%the endend value is endbead-z
- data = reshape(data'(:), nb, 3, []);
+ nsnap = size(data, 1)/nb;
+ data = permute(reshape(data(:,:)', 3, nb, nsnap), [3, 2, 1]);
+
+ % rdata structure is
+ % rdata ( isnaphot, ibead, dim )
+
% the number of snapshots
- nsnap = size(data, 3);
- extfun = zeros(nsnap, 1);
-%the end-to-end distance
-extfun=sqrt(getend(data,nb,dim,nsnap));
-% extfun(:) = max(data(:, 1, :)) - min(data(:, 1, :));
+ %the end-to-end distance
+ Rhead = squeeze(data(:, 1, :));
+ Rtail = squeeze(data(:, end, :));
+ Re2e = Rtail - Rhead;
+ extfun= sqrt(sumsq (Re2e, 2));
endfunction
-%Get the end-to-end distance over all snapshots
-function [endtoend]=getend(data,nb,dim,nsnap)
-endtoend=0;
-%(x1-x2)^2
-for j=1:nsnap
-for i=1:dim
-endtoend+=(data(i,1,j)-data(nb-3+i,dim,j)).^2;
-end
-end
-endtoend=endtoend/nsnap;
+
+function Rg2 = getrg2(data)
+ nb = size(data, 2);
+ Rcom = mean(data, 2);
+ Rincm = data - repmat(Rcom, [1, nb, 1]);
+ Rg2 = mean(sumsq(Rincm, 3), 2);
endfunction
+
% get the number of beads
function [nbead] = getnbead(fname)
fid = fopen (fname);
Oops, something went wrong.

0 comments on commit 7b529dc

Please sign in to comment.