6
variable type string "_charge"
10
read_data cnt_9_0_100${type}.init
11
lattice diamond 3.6 # NOTE ???
13
pair_coeff * * ./CH.airebo C
16
# PARAMETERS-----------------------------
17
variable dt equal 0.0005
18
variable L equal zhi-zlo
19
variable zhi equal zhi
20
variable zTip equal ${zhi}-2.0 # 2 4
21
variable zFree equal zhi
22
variable R equal 12.1/2
23
variable xhiFE equal 5.0*$R
24
variable xloFE equal -${xhiFE}
25
variable yhiFE equal $R
26
variable yloFE equal -${yhiFE}
27
variable zloFE equal zlo+10 # create fixed ghosts
28
variable zhiFE equal zhi+(zhi-${zloFE})/12*2
29
variable Lfree equal zhi-${zloFE}
30
variable nx equal 10 # 5
31
variable nz equal 14 # 12
32
print "Length $L [${zloFE}, ${zhiFE}] ${zTip}"
34
#variable E equal 0.1 # bias 1.0
35
variable Vb equal 0.1 # 0.1 #0.5 #0.0 # bias
36
variable Vg equal 0.5 # 1.0 # 5.0 0.5 #50.0 # 0.5 # gate/modulation
37
print "bias voltage ${Vb}, gate voltage ${Vg}"
39
variable ng equal 20 # 80 # 10
40
variable nb equal 2 # 3
42
variable n equal 100000
44
# END -----------------------------------
46
region TIP block INF INF INF INF ${zTip} INF units box
49
#region FIXED block INF INF INF INF INF ${zLoFE} units box
50
#group FIXED region FIXED
52
#group FREE subtract all FIXED
54
region feRegion block ${xloFE} ${xhiFE} ${yloFE} ${yhiFE} ${zloFE} ${zhiFE} units box
55
group internal region feRegion
56
group FIXED subtract all internal
57
fix FIX FIXED setforce 0 0 0
59
variable nAll equal count(all)
60
variable nGhost equal count(all)-count(internal)
61
variable nTip equal count(TIP)
62
print ">>> number of stationary ghosts: ${nGhost} of ${nAll}"
63
print ">>> number of tip atoms : ${nTip}"
66
neigh_modify every 10 delay 0 check no
69
variable tag string "electrostatic_bending"
72
variable C equal -0.025
74
variable c equal $C/${nTip}
75
set group TIP charge $c
78
fix AtC internal atc electrostatic CNT_id.mat
79
fix_modify AtC include atomic_charge
80
fix_modify AtC internal_quadrature off
81
# note weights don't affect phi or f
82
fix_modify AtC atom_weight constant internal 1.0
83
fix_modify AtC extrinsic short_range off
84
#fix_modify AtC atom_element_map eulerian 1
85
fix_modify AtC control momentum flux
87
fix_modify AtC mesh create ${nx} 1 ${nz} feRegion f p f
90
fix_modify AtC initial displacement x all 0.0
91
fix_modify AtC initial displacement y all 0.0
92
fix_modify AtC initial displacement z all 0.0
93
fix_modify AtC initial velocity x all 0.0
94
fix_modify AtC initial velocity y all 0.0
95
fix_modify AtC initial velocity z all 0.0
96
fix_modify AtC initial electric_potential all 0.0
99
variable t equal 1.1*$R
100
fix_modify AtC mesh create_nodeset tube -$t $t -$t $t ${zloFE} ${zFree} units box
101
fix_modify AtC mesh create_nodeset lefttube -$t $t -$t $t ${zloFE} ${zloFE} units box
102
fix_modify AtC mesh create_nodeset rbc INF INF INF INF ${zhiFE} ${zhiFE} units box
103
fix_modify AtC mesh create_nodeset lbc INF INF INF INF ${zloFE} ${zloFE} units box
104
fix_modify AtC mesh create_nodeset top ${xhiFE} ${xhiFE} INF INF INF INF units box
105
fix_modify AtC mesh create_nodeset bot ${xloFE} ${xloFE} INF INF INF INF units box
107
# boundary conditions
108
fix_modify AtC fix displacement x lbc 0.
109
fix_modify AtC fix displacement y lbc 0.
110
fix_modify AtC fix displacement z lbc 0.
111
fix_modify AtC fix velocity x lbc 0.
112
fix_modify AtC fix velocity y lbc 0.
113
fix_modify AtC fix velocity z lbc 0.
115
fix_modify AtC fix electric_potential lbc 0
117
fix_modify AtC fix electric_potential rbc ${Vb}
119
fix_modify AtC fix electric_potential bot ${Vg}
123
compute q all property/atom q
124
compute Q all reduce sum c_q
125
compute FSUM all reduce sum fx fy fz
126
compute RSUM internal reduce sum fx fy fz
128
thermo_style custom step etotal ke c_CM[1] c_CM[2] c_CM[3] &
129
c_Q f_AtC[4] f_AtC[5] f_AtC[6] f_AtC[7] f_FIX[1] f_FIX[2] f_FIX[3] f_AtC c_FSUM[1] c_RSUM[1]
131
fix_modify AtC output ${tag}FE 100000000 full_text # $s full_text #binary
132
fix_modify AtC output index step
133
# NOTE not recognized as vector by paraview
134
variable uX atom x-f_AtC[1]
135
variable uY atom y-f_AtC[2]
136
variable uZ atom z-f_AtC[3]
137
variable rho atom mass*f_AtC[4]
138
dump CONFIG all custom $s ${tag}.dmp id type x y z v_uX v_uY v_uZ v_rho
143
# [eV/A * A^2] --> [N m]
144
variable eV2J equal 1.60217646e-19
145
variable A2m equal 1.e-10
148
min_modify line quadratic
149
variable Vg equal 0.1
150
variable Lx equal 1.0
152
#compute RSUM FREE reduce sum fx fy fz
153
#dump CONFIG all custom 10000 ${tag}.dmp id type x y z c_U[1] c_U[2] c_U[3] fx fy fz
155
variable i loop ${ng}
157
variable b equal ($i-1)*${Vg}/(${ng}-1)/${Lx}
159
fix_modify AtC fix electric_potential all linear 0 0 0 $b 0 $a 0
161
min_modify line quadratic
162
#minimize 0 0 100000 100000
163
minimize 0 0 1000 1000
165
min_modify line backtrack
166
#minimize 0 0 100000 100000
167
minimize 0 0 1000 1000
168
fix_modify AtC output now
170
# u = F L^3 / 3 EI --> EI = F L^3 / 3 u
171
variable u equal c_CM[1]
172
variable uz equal c_CM[3]
173
# variable F equal f_AtC[5]
174
# variable Fz equal f_AtC[7]
175
variable F equal c_RSUM[1]
176
variable Fz equal c_RSUM[3]
177
variable R equal $F-$C*$b
178
variable Rz equal ${Fz}-$C*$a
179
variable EI equal $F*${Lfree}*${Lfree}*${Lfree}/3./$u
180
variable EI equal ${EI}*${eV2J}*${A2m}
181
#print "flexural rigidity ${EI} [Nm^2] NOTE z force"
183
print ">> V $b $a F $F ${Fz} u $u ${uz} c $c phi 0 EI ${EI} R $R ${Rz}"