echo both
units		metal
atom_style charge
dielectric      1.

boundary	s s f
# read in CNT
read_data  min_CNT_dos.data
set group all charge 0
lattice       diamond  3.6 
pair_style  airebo 3.0
pair_coeff  * * ./CH.airebo C
mass          *  12.01

compute q all property/atom q
compute Q  all reduce sum c_q

# PARAMETERS-----------------------------
# [eV/A * A^2] --> [N m]
variable eV2J equal 1.60217646e-19
variable A2m  equal 1.e-10

variable Lx equal xhi-xlo
variable L equal zhi-zlo
variable zTip equal zhi-3.5 
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 Vb equal 0.1 
variable Vg equal 0.15 
variable V0 equal 1. # 2.
print "bias voltage ${Vb}, gate voltage ${Vg}"

variable ng equal 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	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

thermo 10
set group all image 0 0 0
compute CM TIP com
thermo_style    custom step c_Q etotal c_CM[1] c_CM[3]
#minimize 0 0 1000 1000
#write_restart min_CNT_dos.rst
run 0
#EXIT
variable L equal c_CM[1]
variable Lx  equal $L
variable dx  equal c_CM[1]-${Lx}
variable L equal c_CM[3]
variable Lz  equal $L
variable dz  equal c_CM[3]-${Lz}
print "initial ${Lx} ${Lz} "

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

neighbor	5. bin
neigh_modify	every 10 delay 0 check no
timestep        0
thermo 100

# coupling ............................................................
fix         AtC internal   atc electrostatic CNT_electrostatic2.mat
fix_modify  AtC omit atomic_charge
#fix_modify  AtC internal_quadrature off ## NOTE active -> error 
# note weights don't affect phi or f i.e. they divide out
fix_modify  AtC atom_weight constant internal 1.0
fix_modify  AtC extrinsic short_range off
fix_modify  AtC source_integration atom
fix_modify  AtC atom_element_map eulerian 1
fix_modify  AtC control momentum none # flux

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

# 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.

# minimize .............................................................
compute FSUM all reduce sum fx fy fz
compute RSUM internal reduce sum fx fy fz

thermo          $s
fix_modify      AtC  output  electrostatic_bending_dosFE 100000000 full_text binary
fix_modify      AtC  output  index step      

# store original (reference) coordinates
fix X all store/state 0 x y z

# NOTE not recognized as vector by paraview - due to dump2ensight
variable uX atom x-f_X[1]
variable uY atom y-f_X[2]
variable uZ atom z-f_X[3]
#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 100000 electrostatic_bending_dos.dmp id type x y z c_q v_uX v_uY v_uZ v_rho

reset_timestep 0
log electrostatic_bending_dos.log

thermo 10 # 1 # 10
min_modify      line quadratic
variable a equal 0
variable i loop ${ng}
thermo_style    custom step c_Q pe v_dx v_dz f_FIX[1] f_FIX[3]
label loop_i
  variable b equal ($i-1)*${Vg}/${ng}/${Lz}
  fix_modify AtC  fix electric_potential all linear 0 0 0 $b 0 $a ${V0} # <<<ALL
  min_style cg
  min_modify line quadratic 
  minimize 0 0 1000 1000
  #min_style sd
  #min_modify line backtrack 
  #minimize 0 0 1000 1000
  fix_modify AtC output now 

  # u = F L^3 / 3 EI --> EI = F L^3 / 3 u
  variable Q equal c_Q
  variable ux equal ${dx}
  variable uz equal ${dz}
  variable Fx equal f_FIX[1]
  variable Fz equal f_FIX[3]
  variable EI equal ${Fx}*${Lfree}*${Lfree}*${Lfree}/3./${ux}
  variable EI equal ${EI}*${eV2J}*${A2m}
  print ">> V $b $a F ${Fx} ${Fz} u ${ux} ${uz} Q $Q EI ${EI}"

  next i
jump SELF loop_i