-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgcmc.in
89 lines (69 loc) · 3.02 KB
/
gcmc.in
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
units lj
atom_style full
boundary p p f
processors * * 1
### Force field ###
# r_cut
pair_style lj/cut 2.5
pair_modify shift yes
bond_style fene
special_bonds fene
read_data ${initial_data_file}
# K R0 e s
bond_coeff 1 30.0 1.5 1.0 1.0
# 'Static' pair coefficients
# e s cutoff
pair_coeff * * 1.0 1.0 1.12446 # General repulsive
pair_coeff 3 3 1.0 1.0 2.5 # Solvent-solvent
# Dynamic pair coefficents
# e s cutoff
pair_coeff 1*2 1*2 ${epp} 1.0 2.5
pair_coeff 1*2 3 ${eps} 1.0 ${cps}
### Run settings ###
timestep 0.015
run_style respa 2 2
neighbor 0.3 bin
neigh_modify delay 1 every 1 check yes
# Grafting layer
group grafting type 1
group mobile subtract all grafting
fix walls all wall/harmonic zlo EDGE 100 1 1 zhi EDGE 100 1 1
fix freezeGrafting grafting setforce 0 0 0
neigh_modify exclude group grafting grafting
# GCMC
region solvent block EDGE EDGE EDGE EDGE $(zhi-v_sol_height-1) $(zhi-1) units box
group solvent type 3
group polymer type 2
variable Nvap equal count(solvent)
fix mc solvent gcmc 10000 1000 0 3 12843 ${T_RT} ${mu} 0 region solvent group solvent pressure $p
compute vaptemp solvent temp
compute_modify vaptemp dynamic/dof yes
fix vapnvt solvent nvt temp ${T_RT} ${T_RT} $(10*dt)
fix_modify vapnvt temp vaptemp
fix polynvt polymer nvt temp ${T_RT} ${T_RT} $(10*dt)
### Output ###
# Density profiles
compute cc all chunk/atom bin/1d z 0.0 0.25
fix prof_dens_poly polymer ave/chunk 100 100 10000 cc density/number file PolyDens.dat
fix prof_dens_solv solvent ave/chunk 100 100 10000 cc density/number file SolvDens.dat
# Temperature profile
fix prof_temp all ave/chunk 100 100 10000 cc temp file Temp.dat
# Pressure profile
compute atom_stress all stress/atom thermo_temp
variable atom_press atom -(c_atom_stress[1]+c_atom_stress[2]+c_atom_stress[3])/3
variable atom_press_z atom -c_atom_stress[3]
fix prof_press solvent ave/chunk 100 100 10000 cc v_atom_press file Pressure.dat
fix prof_press_z solvent ave/chunk 100 100 10000 cc v_atom_press_z file PressureZ.dat
# Output pressure in solvent region
variable solvol equal lx*ly*v_sol_height
compute peratom solvent stress/atom vaptemp
compute p solvent reduce/region solvent sum c_peratom[1] c_peratom[2] c_peratom[3]
variable press equal -(c_p[1]+c_p[2]+c_p[3])/(3*v_solvol)
fix sol_press solvent ave/time 100 100 10000 v_press file SolPress.dat
compute Tmobile mobile temp
thermo 10000
thermo_style custom step temp pe etotal press v_press c_Tmobile v_Nvap cpuremain
shell mkdir img_gcmc
dump 1 all image 100000 img_gcmc/*.jpg type type view 90 90 size 1000 1000 zoom 2
run $r
write_data gcmc.data