echo both
units		metal
atom_style charge
dielectric      1.

variable type string "_charge"

boundary	s s f
# read in CNT
read_data  cnt_9_0_100${type}.data
lattice       diamond  3.6 # NOTE ???
pair_style  airebo 3.0
pair_coeff  * * ./CH.airebo C
mass          *  12.01

# PARAMETERS-----------------------------
variable dt equal 0.0005
variable L equal zhi-zlo
variable zhi equal zhi
variable zTip equal ${zhi}-2.0 # 2 4
variable zFree equal zhi
variable R equal 12.1/2
variable xhiFE equal 5.0*$R
variable xloFE equal -${xhiFE}
variable yhiFE equal $R
variable yloFE equal -${yhiFE}
variable zloFE equal zlo+10 # create fixed ghosts
variable zhiFE equal zhi+(zhi-${zloFE})/12*2
variable Lfree equal zhi-${zloFE}
variable nx equal 10 # 5
variable nz equal 14 # 12
print "Length $L [${zloFE}, ${zhiFE}] ${zTip}"

#variable E  equal 0.1  # bias 1.0
variable Vb equal 0.1 # 0.1 #0.5 #0.0 # bias
variable Vg equal 0.5 # 1.0 # 5.0 0.5 #50.0 # 0.5 # gate/modulation
print "bias voltage ${Vb}, gate voltage ${Vg}"

variable ng equal 20 # 80 # 10
variable nb equal 2 # 3

variable n equal 100000
variable s equal 250
# END -----------------------------------

region	TIP block INF INF INF INF ${zTip} INF  units box
group	TIP region TIP

#region  FIXED block INF INF INF INF INF ${zLoFE}  units box
#group   FIXED region FIXED

#group   FREE subtract all FIXED

region	feRegion block ${xloFE} ${xhiFE} ${yloFE} ${yhiFE} ${zloFE} ${zhiFE} units box
group internal region feRegion
group FIXED subtract all internal
fix FIX FIXED setforce 0 0 0

variable nAll   equal count(all)
variable nGhost equal count(all)-count(internal)
variable nTip   equal count(TIP)
print ">>> number of stationary ghosts:  ${nGhost} of ${nAll}"
print ">>> number of tip atoms :  ${nTip}"

neighbor	5. bin
neigh_modify	every 10 delay 0 check no
timestep        ${dt} 
thermo 100
variable tag string "electrostatic_bending"

# set charge on tip
variable C equal -0.025 
print "charge $C [e]"
variable c equal $C/${nTip}
set group TIP charge $c

# coupling
fix         AtC internal   atc electrostatic CNT_id.mat
fix_modify  AtC include atomic_charge
fix_modify  AtC internal_quadrature off
# note weights don't affect phi or f
fix_modify  AtC atom_weight constant internal 1.0
fix_modify  AtC extrinsic short_range off
#fix_modify  AtC atom_element_map eulerian 1
fix_modify  AtC control  momentum flux

fix_modify  AtC mesh create ${nx} 1 ${nz} feRegion f p f

# initial conditions
fix_modify  AtC  initial displacement x all 0.0
fix_modify  AtC  initial displacement y all 0.0
fix_modify  AtC  initial displacement z all 0.0
fix_modify  AtC  initial velocity x all 0.0
fix_modify  AtC  initial velocity y all 0.0
fix_modify  AtC  initial velocity z all 0.0
fix_modify  AtC  initial electric_potential all 0.0

# node sets
variable t equal 1.1*$R
fix_modify  AtC mesh create_nodeset tube -$t $t -$t $t ${zloFE} ${zFree} units box
fix_modify  AtC mesh create_nodeset lefttube -$t $t -$t $t ${zloFE} ${zloFE} units box
fix_modify  AtC mesh create_nodeset rbc INF INF INF INF ${zhiFE} ${zhiFE} units box
fix_modify  AtC mesh create_nodeset lbc INF INF INF INF ${zloFE} ${zloFE} units box
fix_modify  AtC mesh create_nodeset top ${xhiFE} ${xhiFE} INF INF INF INF units box
fix_modify  AtC mesh create_nodeset bot ${xloFE} ${xloFE} INF INF INF INF units box

# boundary conditions
fix_modify AtC  fix displacement x lbc 0.
fix_modify AtC  fix displacement y lbc 0.
fix_modify AtC  fix displacement z lbc 0.
fix_modify AtC  fix velocity     x lbc 0.
fix_modify AtC  fix velocity     y lbc 0.
fix_modify AtC  fix velocity     z lbc 0.
# ground
fix_modify AtC  fix electric_potential lbc 0
# bias
fix_modify AtC  fix electric_potential rbc ${Vb}
# gate
fix_modify AtC  fix electric_potential bot ${Vg} 

# run
compute CM TIP com
compute q all property/atom q
compute Q  all reduce sum c_q
compute FSUM all reduce sum fx fy fz
compute RSUM internal reduce sum fx fy fz

thermo_style    custom step etotal ke c_CM[1] c_CM[2] c_CM[3]  &
  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]
thermo          $s
fix_modify      AtC  output  ${tag}FE 100000000 full_text # $s full_text #binary
fix_modify      AtC  output  index step      
# NOTE not recognized as vector by paraview
variable uX atom x-f_AtC[1]
variable uY atom y-f_AtC[2]
variable uZ atom z-f_AtC[3]
variable rho atom mass*f_AtC[4]
dump  CONFIG all custom $s ${tag}.dmp id type x y z v_uX v_uY v_uZ v_rho

reset_timestep 0
log ${tag}.log

# [eV/A * A^2] --> [N m]
variable eV2J equal 1.60217646e-19
variable A2m  equal 1.e-10
thermo 10
timestep   0.0
min_modify      line quadratic
variable Vg equal 0.1
variable Lx equal 1.0
variable ng equal 4
#compute RSUM FREE reduce sum fx fy fz
#dump  CONFIG all custom 10000 ${tag}.dmp id type x y z c_U[1] c_U[2] c_U[3] fx fy fz
variable a equal 0
  variable i loop ${ng}
  label loop_i
    variable b equal ($i-1)*${Vg}/(${ng}-1)/${Lx}

    fix_modify AtC  fix electric_potential all linear 0 0 0 $b 0 $a 0
    min_style cg
    min_modify line quadratic 
    #minimize 0 0 100000 100000
    minimize 0 0 1000 1000
    min_style sd
    min_modify line backtrack 
    #minimize 0 0 100000 100000
    minimize 0 0 1000 1000
    fix_modify AtC output now 

    # u = F L^3 / 3 EI --> EI = F L^3 / 3 u
    variable u equal c_CM[1]
    variable uz equal c_CM[3]
#   variable F equal f_AtC[5]
#   variable Fz equal f_AtC[7]
    variable F equal c_RSUM[1]
    variable Fz equal c_RSUM[3]
    variable R equal $F-$C*$b
    variable Rz equal ${Fz}-$C*$a
    variable EI equal $F*${Lfree}*${Lfree}*${Lfree}/3./$u
    variable EI equal ${EI}*${eV2J}*${A2m}
    #print "flexural rigidity ${EI}  [Nm^2]  NOTE z force"

    print ">> V $b $a F $F ${Fz} u $u ${uz} c $c phi 0 EI ${EI} R $R ${Rz}"
    next i
  jump SELF loop_i