diff --git a/examples/USER/atc/elastic/in.cnt_electrostatic b/examples/USER/atc/elastic/in.cnt_electrostatic
deleted file mode 100644
index ee4847c1a..000000000
--- a/examples/USER/atc/elastic/in.cnt_electrostatic
+++ /dev/null
@@ -1,107 +0,0 @@
-# needs description
-#  E = - grad \phi
-#  f = q E 
-echo both
-units		metal
-atom_style	atomic
-
-lattice       diamond  3.6
-pair_style    tersoff
-boundary	s s f
-
-read_data  tube_8_4.init
-
-# PARAMETERS-----------------------------
-variable L equal zhi-zlo
-variable R equal 12.1/2
-variable xhiFE equal 5.0*$R
-variable xloFE equal -${xhiFE}
-variable yhiFE equal $R
-variable yloFE equal -${xhiFE}
-variable zhiFE equal zhi
-variable zloFE equal zlo+10
-print "Length $L [${zloFE}, ${zhiFE}]"
-
-variable E equal 10.0
-print "Electric field $E"
-variable drhodx equal 0.0001
-variable s equal 50
-# END -----------------------------------
-
-#pair_coeff    * * SiC.tersoff  C
-pair_coeff    * * ../../../../potentials/SiC.tersoff  C
-mass          *  12.01
-
-# all atoms simulation
-region		feRegion block ${xloFE} ${xhiFE} ${yloFE} ${yhiFE} ${zloFE} ${zhiFE} units box
-group		internal region feRegion
-
-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.0005
-
-# coupling
-fix         AtC internal   atc electrostatic CNT_electrostatic.mat
-fix_modify  AtC  internal_quadrature off
-fix_modify	AtC  omit atomic_charge 
-fix_modify  AtC mesh create 5  1  12  feRegion f p f
-
-# initial & boundary 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
-
-variable a equal -$R-0.1
-variable b equal  $R+0.1
-fix_modify  AtC mesh create_nodeset tube $a $b $a $b ${zloFE} ${zhiFE} units box
-variable a equal ${zloFE}-0.1
-variable b equal ${zloFE}+0.1
-fix_modify  AtC mesh create_nodeset lbc ${xloFE} ${xhiFE} ${xloFE} ${xhiFE} $a $b units box
-variable a equal ${xhiFE}-0.1
-variable b equal ${xhiFE}+0.1
-fix_modify  AtC mesh create_nodeset top $a $b ${yloFE} ${yhiFE} ${zloFE} ${zhiFE} units box
-variable a equal ${xloFE}-0.1
-variable b equal ${xloFE}+0.1
-fix_modify  AtC mesh create_nodeset bot $a $b ${yloFE} ${yhiFE} ${zloFE} ${zhiFE} units box
-
-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.
-fix_modify AtC  fix electron_density all 0.0
-#fix_modify  AtC  fix electron_density tube 0.2
-fix_modify AtC  fix electron_density tube linear 0 0 0 0 0 ${drhodx} 0  
-fix_modify AtC  fix electric_potential all linear 0 0 0 $E 0 0 0 0
-
-fix_modify      AtC control  momentum flux
-
-# run
-thermo_style    custom step cpu etotal ke
-thermo          1#$s
-fix_modify      AtC  output        cnt_electrostaticFE $s full_text
-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 cnt_electrostatic.dmp id type x y z v_uX v_uY v_uZ v_rho
-log cnt_electrostatic.log
-#run   	        1000
-
-# fixed charge, bc on potential
-fix_modify AtC  unfix electric_potential all
-fix_modify AtC  fix electric_potential lbc linear 0 0 0 $E 0 0 0 0
-fix_modify AtC  fix electric_potential top linear 0 0 0 $E 0 0 0 0
-fix_modify AtC  fix electric_potential bot linear 0 0 0 $E 0 0 0 0
-run   	        1000
diff --git a/examples/USER/atc/elastic/in.cnt_electrostatic2 b/examples/USER/atc/elastic/in.cnt_electrostatic2
deleted file mode 100644
index 7edd8e0f7..000000000
--- a/examples/USER/atc/elastic/in.cnt_electrostatic2
+++ /dev/null
@@ -1,118 +0,0 @@
-# needs description
-#  E = - grad \phi
-#  f = q E 
-
-# NOTE tangent is constant for LAGRANGIAN but not exact...
-# NOTE try one atom and one free node
-
-# issue is one of magnitude of E since
-# tangent = perm BB - dn/dphi NN
-# CHECK CONDITIONING?
-
-echo both
-units		metal
-atom_style	atomic
-
-lattice       diamond  3.6
-boundary      f f f
-
-#region  box block -4.3458792459312239 4.3458792459312328 -4.3458792459312168 4.3458792459312203 0.0 104.29551668 units box
-#region  box block 0 10 0 10 0 10
-#create_box 1 box
-#group   box region box
-#atom_modify  sort 0 1
-
-pair_style    tersoff
-read_data  tube_8_4.data
-pair_coeff    * * ../../../../potentials/SiC.tersoff  C
-mass          *  12.01
-
-# PARAMETERS-----------------------------
-variable L equal zhi-zlo
-variable R equal 12.1/2
-variable xhiFE equal 5.0*$R
-variable xloFE equal -${xhiFE}
-variable yhiFE equal $R
-variable yloFE equal -${xhiFE}
-variable zhiFE equal zhi
-variable zloFE equal zlo+10
-print "Length $L [${zloFE}, ${zhiFE}]"
-
-variable E equal  0.01  # 1.0 10.0 0.01
-variable V equal  $E*${zloFE} 
-variable V equal  2
-print "Electric field $E ref.voltage $V"
-variable s equal  20
-# END -----------------------------------
-
-# all atoms simulation
-region		feRegion block ${xloFE} ${xhiFE} ${yloFE} ${yhiFE} ${zloFE} ${zhiFE} units box
-group		internal region feRegion
-
-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.0005
-
-# coupling
-fix         AtC internal   atc electrostatic-equilibrium CNT_electrostatic2.mat
-fix_modify  AtC omit atomic_charge 
-#fix_modify  AtC internal_quadrature off <<<< ???
-variable alat equal 1.42
-variable w equal ${alat}*${alat}*3.*sqrt(3.)/4.0
-#variable w equal 10
-fix_modify  AtC atom_weight constant internal $w
-fix_modify  AtC source_integration atom
-fix_modify  AtC atom_element_map eulerian 1
-fix_modify  AtC mesh create 5  1  12  feRegion f p f
-fix_modify  AtC control  momentum flux
-#fix_modify  AtC extrinsic poisson_solver max_iterations 5
-
-# 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
-
-#variable a equal -$R-0.1
-#variable b equal  $R+0.1
-#fix_modify  AtC mesh create_nodeset tube $a $b $a $b ${zloFE} ${zhiFE} units box
-fix_modify  AtC mesh create_nodeset lbc ${xloFE} ${xhiFE} ${xloFE} ${xhiFE} ${zloFE} ${zloFE} units box
-fix_modify  AtC mesh create_nodeset rbc ${xloFE} ${xhiFE} ${xloFE} ${xhiFE} ${zhiFE} ${zhiFE} units box
-fix_modify  AtC mesh create_nodeset bot ${xloFE} ${xloFE} ${yloFE} ${yhiFE} ${zloFE} ${zhiFE} units box
-fix_modify  AtC mesh create_nodeset top ${xhiFE} ${xhiFE} ${yloFE} ${yhiFE} ${zloFE} ${zhiFE} 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.
-fix_modify AtC  fix electric_potential lbc linear 0 0 0 $E 0 0 $V
-#fix_modify AtC  fix electric_potential rbc linear 0 0 0 $E 0 0 $V
-##fix_modify AtC  fix electric_potential top linear 0 0 0 $E 0 0 $V
-fix_modify AtC  fix electric_potential bot linear 0 0 0 $E 0 0 $V
-
-# run
-thermo_style    custom step cpu etotal ke
-thermo          $s
-fix_modify      AtC  output cnt_electrostatic2FE $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 cnt_electrostatic2.dmp id type x y z v_uX v_uY v_uZ v_rho
-log cnt_electrostatic2.log
-
-#run   	        $s
-run   	        100
-
-# NOTE try fix charge on tip
diff --git a/examples/USER/atc/elastic/in.cnt_fixed_charge b/examples/USER/atc/elastic/in.cnt_fixed_charge
deleted file mode 100644
index 99e2d54ad..000000000
--- a/examples/USER/atc/elastic/in.cnt_fixed_charge
+++ /dev/null
@@ -1,145 +0,0 @@
-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}.init
-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
-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.0 # bias
-variable Vg equal 0.5 # gate/modulation
-print "bias voltage ${Vb}, gate voltage ${Vg}"
-
-variable n equal 100000
-variable s equal 250
-# END -----------------------------------
-
-# all atoms simulation
-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 100
-#minimize 0 0 1000 1000
-#write_restart cnt_in_box0.rst
-
-region		TIP block INF INF INF INF ${zTip} INF  units box
-group		TIP region TIP
-
-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} 
-variable tag string "cnt_fixed_charge"
-
-# set charge on tip
-variable C equal -0.1 # -0.01 -0.0001102
-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
-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]
-thermo          $s
-log ${tag}.log
-#run   $n 
-#run   $n
-thermo 10
-timestep   0.0
-min_modify      line quadratic
-minimize 0 0 1000 1000
-
-# u = F L^3 / 3 EI --> EI = F L^3 / 3 u
-variable u equal c_CM[1]
-variable F equal f_AtC[5]
-# [eV/A * A^2] --> [N m]
-variable eV2J equal 1.60217646e-19
-variable A2m  equal 1.e-10
-variable EI equal $F*${Lfree}*${Lfree}*${Lfree}/3./$u
-variable EI equal ${EI}*${eV2J}*${A2m}
-print "flexural rigidity ${EI}  [Nm^2]  NOTE z force"
-# flexural rigidity 6.716732985e-25  [Nm^2]
-
-fix_modify      AtC  output  ${tag}FE 1 full_text 
-fix_modify      AtC  output  index step      
-run 1
diff --git a/examples/USER/atc/elastic/in.electrostatic_bending b/examples/USER/atc/elastic/in.electrostatic_bending
deleted file mode 100644
index da4f03fc7..000000000
--- a/examples/USER/atc/elastic/in.electrostatic_bending
+++ /dev/null
@@ -1,186 +0,0 @@
-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}.init
-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
-
diff --git a/examples/USER/atc/elastic/in.electrostatic_bending_dos b/examples/USER/atc/elastic/in.electrostatic_bending_dos
deleted file mode 100644
index 1a37b6a70..000000000
--- a/examples/USER/atc/elastic/in.electrostatic_bending_dos
+++ /dev/null
@@ -1,164 +0,0 @@
-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
diff --git a/lib/atc/ATC_Coupling.cpp b/lib/atc/ATC_Coupling.cpp
index c5df4d741..b60063883 100644
--- a/lib/atc/ATC_Coupling.cpp
+++ b/lib/atc/ATC_Coupling.cpp
@@ -1,2215 +1,2160 @@
 // ATC headers
 #include "ATC_Coupling.h"
 #include "FE_Engine.h"
 #include "Array.h"
 #include "Array2D.h"
 #include "ATC_Error.h"
 #include "PrescribedDataManager.h"
 #include "AtomicRegulator.h"
 #include "TimeIntegrator.h"
 #include "PhysicsModel.h"
 #include "AtomToMoleculeTransfer.h"
 #include "MoleculeSet.h"
 #include "FieldManager.h"
 
 using std::string;
 using std::map;
 using std::pair;
 using std::set;
 
 namespace ATC {
   //--------------------------------------------------
   ATC_Coupling::ATC_Coupling(string groupName, double ** & perAtomArray, LAMMPS_NS::Fix * thisFix) :
     ATC_Method(groupName, perAtomArray, thisFix),
     consistentInitialization_(false),
     equilibriumStart_(false),
     useFeMdMassMatrix_(false),
     trackCharge_(false),
     temperatureDef_(NONE),
     prescribedDataMgr_(NULL),
     physicsModel_(NULL),
     extrinsicModelManager_(this),
     atomicRegulator_(NULL),
     atomQuadForInternal_(true),
     elementMask_(NULL),
     elementMaskMass_(NULL),
     elementMaskMassMd_(NULL),
     nodalAtomicMass_(NULL),
     nodalAtomicCount_(NULL),
     nodalAtomicHeatCapacity_(NULL),
     internalToMask_(NULL),
     internalElement_(NULL),
     ghostElement_(NULL),
     nodalGeometryType_(NULL),
     bndyIntType_(NO_QUADRATURE),
     bndyFaceSet_(NULL),
     atomicWeightsMask_(NULL),
     shpFcnMask_(NULL),
     shpFcnDerivsMask_(NULL),
     sourceIntegration_(FULL_DOMAIN)
   {
     // size the field mask
     fieldMask_.reset(NUM_FIELDS,NUM_FLUX); 
     fieldMask_ = false;
     // default: no consistent mass matrices
     useConsistentMassMatrix_.reset(NUM_FIELDS);
     useConsistentMassMatrix_ = false;
     mdMassNormalization_ = true;
     // check to see if lammps has any charges
     if (lammpsInterface_->atom_charge()) trackCharge_ = true;
     // default: perform velocity verlet
     integrateInternalAtoms_ = true;
   }
   //--------------------------------------------------
   ATC_Coupling::~ATC_Coupling()
   {
     interscaleManager_.clear(); 
     if (feEngine_) { delete feEngine_; feEngine_ = NULL; } 
     if (physicsModel_) delete physicsModel_;
     if (atomicRegulator_) delete atomicRegulator_;
     if (prescribedDataMgr_) delete prescribedDataMgr_;
     for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
       delete _tiIt_->second;
     }
   }
   //--------------------------------------------------
   // Interactions with LAMMPS fix commands
   // parse input command and pass on to finite element engine 
   //   or physics specific transfers if necessary
   //   revert to physics-specific transfer if no command matches input
   // first keyword is unique to particular class
   // base class keyword matching must apply to ALL physics
   // order:  derived, base, owned objects
   //--------------------------------------------------
   bool ATC_Coupling::modify(int narg, char **arg)
   {
     FieldName thisField;
     int thisIndex;
     int argIdx=0;
 
     bool match = false; 
     
     // gateways to other modules e.g. extrinsic, control, mesh
     // pass off to extrinsic
     if (strcmp(arg[argIdx],"extrinsic")==0) {
       argIdx++;
       match = extrinsicModelManager_.modify(narg-argIdx,&arg[argIdx]);
     }
     // catch special case
     if ((strcmp(arg[argIdx],"control")==0)
       &&(strcmp(arg[argIdx+1],"charge")==0)) {
       match = extrinsicModelManager_.modify(narg-argIdx,&arg[argIdx]);
     }
     // parsing handled here
     else {
       /*! \page man_initial fix_modify AtC initial 
         \section syntax
         fix_modify AtC initial <field> <nodeset> <constant | function>
         - <field> = field name valid for type of physics, temperature | electron_temperature
         - <nodeset> = name of set of nodes to apply initial condition
         - <constant | function> = value or name of function followed by its
           parameters
         \section examples
         <TT> fix_modify atc initial temperature groupNAME 10. </TT>
         \section description
         Sets the initial values for the specified field at the specified nodes.
         \section restrictions
         keyword 'all' reserved in nodeset name
         \section default 
         none
       */
       // set initial conditions
       if (strcmp(arg[argIdx],"initial")==0) {
         argIdx++;
         parse_field(arg,argIdx,thisField,thisIndex);
         string nsetName(arg[argIdx++]);
         XT_Function * f = NULL;
         // parse constant
         if (narg == argIdx+1) {
           f = XT_Function_Mgr::instance()->constant_function(atof(arg[argIdx]));
         }
         // parse function
         else {
           f = XT_Function_Mgr::instance()->function(&(arg[argIdx]),narg-argIdx);
         }
         prescribedDataMgr_->fix_initial_field(nsetName,thisField,thisIndex,f);
         match = true;
       }
 
       /*! \page man_fix_nodes fix_modify AtC fix 
         \section syntax
         fix_modify AtC fix <field> <nodeset> <constant | function>
         - <field> = field name valid for type of physics
         - <nodeset> = name of set of nodes to apply boundary condition
         - <constant | function> = value or name of function followed by its
           parameters
         \section examples
         <TT> fix_modify AtC fix temperature groupNAME 10. </TT> \n
         <TT> fix_modify AtC fix temperature groupNAME 0 0 0 10.0 0 0 1.0 </TT> \n
         \section description
         Creates a constraint on the values of the specified field at specified nodes.
         \section restrictions
         keyword 'all' reserved in nodeset name
         \section related 
         see \ref man_unfix_nodes
         \section default 
         none
       */
       // fix and unfix nodes
       else if (strcmp(arg[argIdx],"fix")==0) {
         argIdx++;
         parse_field(arg,argIdx,thisField,thisIndex);
         string nsetName(arg[argIdx++]);
         XT_Function * f = NULL;
         // fix current value 
         if (narg == argIdx) {
           set<int> nodeSet = (feEngine_->fe_mesh())->nodeset(nsetName);
           set<int>::const_iterator iset;
           const DENS_MAT & field =(fields_.find(thisField)->second).quantity(); 
           for (iset = nodeSet.begin(); iset != nodeSet.end(); iset++) {
             int inode = *iset;
             double v = field(inode,thisIndex);
             f = XT_Function_Mgr::instance()->constant_function(v);
             set<int> one; one.insert(inode);
             prescribedDataMgr_->fix_field(one,thisField,thisIndex,f);
           }
          }
         // parse constant
         else if (narg == argIdx+1) {
           f = XT_Function_Mgr::instance()->constant_function(atof(arg[argIdx]));
           prescribedDataMgr_->fix_field(nsetName,thisField,thisIndex,f);
         }
         // parse function
         else {
           f = XT_Function_Mgr::instance()->function(&(arg[argIdx]),narg-argIdx);
           prescribedDataMgr_->fix_field(nsetName,thisField,thisIndex,f);
         }
         match = true;
       }
 
       /*! \page man_unfix_nodes fix_modify AtC unfix 
         \section syntax
         fix_modify AtC unfix <field> <nodeset> 
         - <field> = field name valid for type of physics
         - <nodeset> = name of set of nodes 
         \section examples
         <TT> fix_modify AtC unfix temperature groupNAME </TT>
         \section description
         Removes constraint on field values for specified nodes.
         \section restrictions
         keyword 'all' reserved in nodeset name
         \section related 
         see \ref man_fix_nodes
         \section default 
         none
       */
       else if (strcmp(arg[argIdx],"unfix")==0) {
         argIdx++;
         parse_field(arg,argIdx,thisField,thisIndex);
         string nsetName(arg[argIdx++]);
         prescribedDataMgr_->unfix_field(nsetName,thisField,thisIndex);
         match = true;
       }
 
     /*! \page man_source fix_modify AtC source
       \section syntax
        fix_modify AtC source <field> <element_set> <value | function>
         - <field> = field name valid for type of physics
         - <element_set> = name of set of elements 
       \section examples
       <TT> fix_modify atc source temperature middle temporal_ramp 10. 0. </TT>
       \section description
       Add domain sources to the mesh. The units are consistent with LAMMPS's 
       units for mass, length and time and are defined by the PDE being solved,
       e.g. for thermal transfer the balance equation is for energy and source
       is energy per time.
       \section restrictions 
       keyword 'all' reserved in element_set name
       \section related
       see \ref man_remove_source
       \section default
       none
     */
       else if (strcmp(arg[argIdx],"source")==0) {
         argIdx++;
         parse_field(arg,argIdx,thisField,thisIndex);
         string esetName(arg[argIdx++]);
         XT_Function * f = NULL;
         // parse constant
         if (narg == argIdx+1) {
           f = XT_Function_Mgr::instance()->constant_function(atof(arg[argIdx]));
         }
         // parse function
         else {
           f = XT_Function_Mgr::instance()->function(&(arg[argIdx]),narg-argIdx);
         }
         prescribedDataMgr_->fix_source(esetName,thisField,thisIndex,f);
         fieldMask_(thisField,PRESCRIBED_SOURCE) = true;
         match = true;
       }
 
     /*! \page man_remove_source fix_modify AtC remove_source
       \section syntax
       fix_modify AtC remove_source <field> <element_set>
         - <field> = field name valid for type of physics
         - <element_set> = name of set of elements 
       \section examples
       <TT> fix_modify atc remove_source temperature groupNAME </TT>
       \section description
       Remove a domain source.
       \section restrictions
       keyword 'all' reserved in element_set name
       \section related
       see \ref man_source
       \section default
     */
       else if (strcmp(arg[argIdx],"remove_source")==0) {
         argIdx++;
         parse_field(arg,argIdx,thisField,thisIndex);
         string esetName(arg[argIdx++]);
         prescribedDataMgr_->unfix_source(esetName,thisField,thisIndex);
         fieldMask_(thisField,PRESCRIBED_SOURCE) = false;
         match = true;
       }
   
     /*! \page man_fix_flux fix_modify AtC fix_flux
       \section syntax
        fix_modify AtC fix_flux <field> <face_set> <value | function>
         - <field> = field name valid for type of physics, temperature | electron_temperature
         - <face_set> = name of set of element faces
       \section examples
        <TT> fix_modify atc fix_flux temperature faceSet 10.0 </TT> \n
        
       \section description
        Command for fixing normal fluxes e.g. heat_flux. 
        This command only prescribes the normal component of the physical flux, e.g. heat (energy) flux.
        The units are in AtC units, i.e. derived from the LAMMPS length, time, and mass scales.
       \section restrictions 
       Only normal fluxes (Neumann data) can be prescribed.
       \section related
       see \ref man_unfix_flux
       \section default
     */
       else if (strcmp(arg[argIdx],"fix_flux")==0) {
         argIdx++;
         parse_field(arg,argIdx,thisField,thisIndex);
         string fsetName(arg[argIdx++]);
         XT_Function * f = NULL;
         // parse constant
         if (narg == argIdx+1) {
           f = XT_Function_Mgr::instance()->constant_function(atof(arg[argIdx]));
         }
         // parse function
         else {
           f = XT_Function_Mgr::instance()->function(&(arg[argIdx]),narg-argIdx);
         }
         prescribedDataMgr_->fix_flux(fsetName,thisField,thisIndex,f);
         fieldMask_(thisField,PRESCRIBED_SOURCE) = true;
         match = true;
       }
 
     /*! \page man_unfix_flux fix_modify AtC unfix_flux
       \section syntax
       fix_modify AtC fix_flux <field> <face_set> <value | function>
         - <field> = field name valid for type of physics, temperature | electron_temperature
         - <face_set> = name of set of element faces
       \section examples
        <TT> fix_modify atc unfix_flux temperature faceSet  </TT> \n
        
       \section description
        Command for removing prescribed normal fluxes e.g. heat_flux, stress. 
       \section restrictions 
       \section related
       see \ref man_unfix_flux
       \section default
     */
       else if (strcmp(arg[argIdx],"unfix_flux")==0) {
         argIdx++;
         parse_field(arg,argIdx,thisField,thisIndex);
         string fsetName(arg[argIdx++]);
         prescribedDataMgr_->unfix_flux(fsetName,thisField,thisIndex);
         fieldMask_(thisField,PRESCRIBED_SOURCE) = false;
         match = true;
       }
 
       
     /*! \page man_fe_md_boundary fix_modify AtC fe_md_boundary 
       \section syntax
       fix_modify AtC fe_md_boundary <faceset | interpolate | no_boundary> [args]
       \section examples
        <TT> fix_modify atc fe_md_boundary interpolate </TT> \n
       \section description
       Specifies different methods for computing fluxes between between the MD and FE integration regions.  Faceset defines a faceset separating the MD and FE regions and uses finite element face quadrature to compute the flux.  Interpolate uses a reconstruction scheme to approximate the flux, which is more robust but less accurate if the MD/FE boundary does correspond to a faceset.  No boundary results in no fluxes between the systems being computed.
       \section restrictions 
       If faceset is used, all the AtC non-boundary atoms must lie within and completely fill the domain enclosed by the faceset.
       \section related
       see \man_boundary_faceset for how to specify the faceset name.
       \section default
       Interpolate.
     */
       else if (strcmp(arg[argIdx],"fe_md_boundary")==0) {
         bndyIntType_ = FE_INTERPOLATION;// default
         if(strcmp(arg[argIdx],"faceset")==0) {
           argIdx++;
           bndyIntType_ = FE_QUADRATURE;
           string name(arg[argIdx++]);
           bndyFaceSet_ = & ( (feEngine_->fe_mesh())->faceset(name));
         } 
         else if (strcmp(arg[argIdx],"interpolate")==0) {
           argIdx++;
           bndyIntType_ = FE_INTERPOLATION;
         }
         else if (strcmp(arg[argIdx],"no_boundary")==0) { 
           bndyIntType_ = NO_QUADRATURE; 
         }
         else {
           throw ATC_Error("Bad boundary integration type");
         }
       }
 
 
 
     /*! \page man_boundary_faceset fix_modify AtC boundary_faceset 
       \section syntax
       fix_modify AtC boundary_faceset <is | add> [args]
       \section examples
       fix_modify AtC boundary_faceset is obndy
       \section description
       This command species the faceset name when using a faceset to compute the MD/FE boundary fluxes.  The faceset must already exist.
       \section restrictions
       This is only valid when fe_md_boundary is set to faceset.
       \section related
       \man_fe_md_boundary
       \section default
     */
       else if (strcmp(arg[argIdx],"boundary_faceset")==0) {
         argIdx++;
         if (strcmp(arg[argIdx],"is")==0) { // replace existing faceset
           argIdx++;
           boundaryFaceNames_.clear();
           string name(arg[argIdx++]);
           boundaryFaceNames_.insert(name);
           match = true;
         }
         else if (strcmp(arg[argIdx],"add")==0) { // add this faceset to list
           argIdx++;
           string name(arg[argIdx]);
           boundaryFaceNames_.insert(name);
           match = true;
         }
       }
 
       /*! \page man_internal_quadrature fix_modify AtC internal_quadrature 
         \section syntax
         fix_modify atc internal_quadrature <on | off> [region]
         \section examples
         <TT> fix_modify atc internal_quadrature off </TT>
         \section description
         Command to use or not use atomic quadrature on internal elements
         fully filled with atoms. By turning the internal quadrature off
         these elements do not contribute to the governing PDE and the fields
         at the internal nodes follow the weighted averages of the atomic data.
         \section optional
         Optional region tag specifies which finite element nodes will be treated
         as being within the MD region.  This option is only valid with
         internal_quadrature off.
         \section restrictions
         \section related 
         \section default 
         on
       */
       else if (strcmp(arg[argIdx],"internal_quadrature")==0) {
         if (initialized_) {
           throw ATC_Error("Cannot change internal_quadrature method after first run");
         }
         argIdx++;
         if (strcmp(arg[argIdx],"on")==0) {
           argIdx++;
           atomQuadForInternal_ = true;
           match = true;
         }
         else if (strcmp(arg[argIdx],"off")==0) {
           argIdx++;
           if (argIdx == narg) {
             atomQuadForInternal_ = false;
             regionID_ = -1;
             match = true;
           }
           else {
             for (regionID_ = 0; regionID_ < lammpsInterface_->nregion(); regionID_++) 
               if (strcmp(arg[argIdx],lammpsInterface_->region_name(regionID_)) == 0) break;
             if (regionID_ < lammpsInterface_->nregion()) {
               atomQuadForInternal_ = false;
               match = true;
             }
             else {
               throw ATC_Error("Region " + string(arg[argIdx]) + " does not exist");
             }
           }
         }
         if (match) {
           needReset_ = true;
         }
       }
 
       else if (strcmp(arg[argIdx],"fix_robin")==0) {
         argIdx++;
         parse_field(arg,argIdx,thisField,thisIndex);
         string fsetName(arg[argIdx++]);
         UXT_Function * f = NULL;
         // parse linear
         if (narg == argIdx+2) {
           f = UXT_Function_Mgr::instance()->linear_function(atof(arg[argIdx]),atof(arg[argIdx+1]));
         }
         // parse function
         else {
           throw ATC_Error("unimplemented function");
         }
         prescribedDataMgr_->fix_robin(fsetName,thisField,thisIndex,f);
         fieldMask_(thisField,ROBIN_SOURCE) = true;
         match = true;
       }
       else if (strcmp(arg[argIdx],"unfix_robin")==0) {
         argIdx++;
         parse_field(arg,argIdx,thisField,thisIndex);
         string fsetName(arg[argIdx++]);
         prescribedDataMgr_->unfix_robin(fsetName,thisField,thisIndex);
         fieldMask_(thisField,ROBIN_SOURCE) = false;
         match = true;
       }
 
       /*! \page man_atomic_charge fix_modify AtC atomic_charge
       \section syntax
       fix_modify AtC <include | omit> atomic_charge
         - <include | omit> = switch to activiate/deactiviate inclusion of intrinsic atomic charge in ATC
       \section examples
        <TT> fix_modify atc compute include atomic_charge </TT>
       \section description
        Determines whether AtC tracks the total charge as a finite element field
       \section restrictions
       Required for:  electrostatics
       \section related
       \section default
       if the atom charge is defined, default is on, otherwise default is off
     */
       else if (strcmp(arg[argIdx],"include")==0) {
         argIdx++;
         if (strcmp(arg[argIdx],"atomic_charge")==0) {
           trackCharge_ = true;
           match = true;
           needReset_ = true;
         }
       }
       else if (strcmp(arg[argIdx],"omit")==0) {
         argIdx++;
         if (strcmp(arg[argIdx],"atomic_charge")==0) {
           trackCharge_ = false;
           match = true;
           needReset_ = true;
         }
       }
 
       /*! \page man_source_integration fix_modify AtC source_integration
       \section syntax
       fix_modify AtC source_integration < fe | atom>
       \section examples
        <TT> fix_modify atc source_integration atom </TT>
       \section description
       \section restrictions
       \section related
       \section default
       Default is fe
     */
       else if (strcmp(arg[argIdx],"source_integration")==0) {
         argIdx++;
         if (strcmp(arg[argIdx],"fe")==0) {
           sourceIntegration_ = FULL_DOMAIN;
         }
-        else {
-          sourceIntegration_ = FULL_DOMAIN_ATOMIC_QUADRATURE_SOURCE;
-        }
         match = true;
       }
 
       /*! \page man_consistent_fe_initialization fix_modify AtC consistent_fe_initialization
       \section syntax
        fix_modify AtC consistent_fe_initialization <on | off>
         - <on|off> = switch to activiate/deactiviate the intial setting of FE intrinsic field to match the projected MD field
       \section examples
        <TT> fix_modify atc consistent_fe_initialization on </TT>
       \section description
       Determines whether AtC initializes FE intrinsic fields (e.g., temperature) to match the projected MD values.  This is particularly useful for fully overlapping simulations.
       \section restrictions
       Can be used with:  thermal, two_temperature.
       Cannot be used with time filtering on. Does not include boundary nodes.
       \section related
       \section default
       Default is off
     */
       else if (strcmp(arg[argIdx],"consistent_fe_initialization")==0) {
         argIdx++;
         if (strcmp(arg[argIdx],"on")==0) {
           if (timeFilterManager_.filter_dynamics())
             throw ATC_Error("Consistent FE initialization cannot be used with time filtering");
           consistentInitialization_ = true;
           match = true;
         }
         else if (strcmp(arg[argIdx],"off")==0) {
           consistentInitialization_ = false;
           match = true;
         }
       }
 
       // switch for equilibrium filtering start
       /*! \page man_equilibrium_start fix_modify AtC equilibrium_start
         \section syntax
         fix_modify AtC equilibrium_start <on|off>
         
         \section examples
         <TT> fix_modify atc equilibrium_start on </TT> \n
         
         \section description
         Starts filtered calculations assuming they start in equilibrium, i.e. perfect finite element force balance.
         
         \section restrictions
         only needed before filtering is begun
         
         \section related
         see \ref man_time_filter
         
         \section default
         on
       */
       else if (strcmp(arg[argIdx],"equilibrium_start")==0) {
         argIdx++;
         if (strcmp(arg[argIdx],"on")==0) {
           equilibriumStart_ = true;
           match = true;
         }
         else if (strcmp(arg[argIdx],"off")==0) {
           equilibriumStart_ = false;
           match = true;
         }
       }
 
       /*! \page man_mass_matrix fix_modify AtC mass_matrix
         \section syntax 
         fix_modify AtC mass_matrix <fe | md_fe>
         - <fe | md_fe> = activiate/deactiviate using the FE mass matrix in the MD region
         \section examples
         <TT> fix_modify atc mass_matrix fe </TT>
         \section description
         Determines whether AtC uses the FE mass matrix based on Gaussian quadrature or based on atomic quadrature in the MD region.  This is useful for fully overlapping simulations to improve efficiency.
         \section restrictions
         Should not be used unless the FE region is contained within the MD region, otherwise the method will be unstable and inaccurate
         \section related
         \section default
         Default is off
       */
 
       else if (strcmp(arg[argIdx],"mass_matrix")==0) {
         argIdx++;
         if (strcmp(arg[argIdx],"fe")==0) {
           useFeMdMassMatrix_ = true;
           match = true;
         }
         else {
           useFeMdMassMatrix_ = false;
           match = true;
         }
         if (match) {
           needReset_ = true;
         }
       }
 
     /*! \page man_material fix_modify AtC material
       \section syntax
       fix_modify AtC material [elementset_name] [material_id]  \n
       \section examples
       <TT> fix_modify AtC material gap_region 2</TT>
       \section description
       Sets the material model in elementset_name to be of type material_id.
       \section restrictions 
       The element set must already be created and the material must be specified in the material file given the the atc fix on construction
       \section related
       \section default
       All elements default to the first material in the material file.
     */
     else if (strcmp(arg[argIdx],"material")==0) {
       argIdx++;
       string elemsetName(arg[argIdx++]);
       int matId = physicsModel_->material_index(arg[argIdx++]);
       using std::set;
       set<int> elemSet = (feEngine_->fe_mesh())->elementset(elemsetName);
       if(elementToMaterialMap_.size() == 0) {
         throw ATC_Error("need mesh before material command");
       }
       // set elementToMaterialMap
       set<int>::const_iterator iset;
       for (iset = elemSet.begin(); iset != elemSet.end(); iset++) {
         int ielem = *iset;
         
         // and the tag a string
         elementToMaterialMap_(ielem) = matId;
       }
       match = true;
       needReset_ = true;
     }
 
     } // end else 
     // no match, call base class parser
     if (!match) {
       match = ATC_Method::modify(narg, arg);
     }
     return match; // return to FixATC
   }
 
   //--------------------------------------------------
    /** PDE type */
    WeakEquation::PDE_Type ATC_Coupling::pde_type(const FieldName fieldName) const
    {
      const WeakEquation * weakEq = physicsModel_->weak_equation(fieldName);
      if (weakEq == NULL) return WeakEquation::PROJECTION_PDE; 
      return weakEq->type();
    }
   //--------------------------------------------------
    /** is dynamic PDE */
    bool ATC_Coupling::is_dynamic(const FieldName fieldName) const
    {
      const WeakEquation * weakEq = physicsModel_->weak_equation(fieldName);
      if (weakEq == NULL) return false; 
      return (physicsModel_->weak_equation(fieldName)->type() == WeakEquation::DYNAMIC_PDE);
    }
 
 
   //--------------------------------------------------
   /** allow FE_Engine to construct data manager after mesh is constructed */
   void ATC_Coupling::construct_prescribed_data_manager (void) {
     prescribedDataMgr_ = new PrescribedDataManager(feEngine_,fieldSizes_);
   }
 
   //--------------------------------------------------
   // pack_fields
   //   bundle all allocated field matrices into a list
   //   for output needs
   //--------------------------------------------------
   void ATC_Coupling::pack_fields(RESTART_LIST & data)
   {
     ATC_Method::pack_fields(data);
     for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
       (_tiIt_->second)->pack_fields(data);
     }
   }
 
   //--------------------------------------------------------------
   // create_physics_model
   // - method to create physics model 
   //--------------------------------------------------------------
   void ATC_Coupling::create_physics_model(const PhysicsType & physicsType,
                                           string matFileName)
   {
     if (physicsModel_) {
       throw ATC_Error("Attempted to create PhysicsModel multiple times in ATC_Coupling");
     }
     // Create PhysicsModel based on physicsType
     switch (physicsType) {
     case NO_PHYSICS :
       break;
     case THERMAL :
       physicsModel_ = new PhysicsModelThermal(matFileName);
       break;
     case ELASTIC :
       physicsModel_ = new PhysicsModelElastic(matFileName);
       break;
     case SHEAR:
       physicsModel_ = new PhysicsModelShear(matFileName);
       break;
     case SPECIES :
       physicsModel_ = new PhysicsModelSpecies(matFileName);
       break;
     case THERMO_ELASTIC :
       physicsModel_ = new PhysicsModelThermoElastic(matFileName);
       break;
     default:
       throw ATC_Error("Unknown physics type in ATC_Coupling::create_physics_model");
     }
   }
 
   //--------------------------------------------------------
   //  construct_methods
   //    have managers instantiate requested algorithms
   //    and methods
   //--------------------------------------------------------
   void ATC_Coupling::construct_methods()
   {
     ATC_Method::construct_methods();
 
     // construct needed time filters for mass matrices
     if (timeFilterManager_.need_reset()) {
       init_filter();
       map<FieldName,int>::const_iterator field;
       for (field = fieldSizes_.begin(); field!=fieldSizes_.end(); field++) {
         FieldName thisField = field->first;
         // fill in mass matrix time filters if needed
         if (!massMatTimeFilters_[thisField])
           massMatTimeFilters_[thisField] = timeFilterManager_.construct(TimeFilterManager::INSTANTANEOUS);
       }
     }
 
     for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
       (_tiIt_->second)->construct_methods();
     }
     atomicRegulator_->construct_methods();
   }
   //-------------------------------------------------------------------
   void ATC_Coupling::init_filter()
   {
     if (timeFilterManager_.need_reset()) {
       map<FieldName,int>::const_iterator field;
       for (field = fieldSizes_.begin(); field!=fieldSizes_.end(); field++) {
         FieldName thisField = field->first;
         int thisSize = field->second;
         (nodalAtomicFieldsRoc_[thisField].set_quantity()).reset(nNodes_,thisSize);
       }
     }
   }
   //--------------------------------------------------------
   void ATC_Coupling::set_fixed_nodes()
   {
     // set fields
     prescribedDataMgr_->set_fixed_fields(time(),
       fields_,dot_fields_,ddot_fields_,dddot_fields_);
 
 
     // set related data
     map<FieldName,int>::const_iterator field;
     for (field = fieldSizes_.begin(); field!=fieldSizes_.end(); field++) {
       FieldName thisField = field->first;
       int thisSize = field->second;
       DENS_MAT & rhs(rhs_[thisField].set_quantity());
       for (int inode = 0; inode < nNodes_ ; ++inode) {
         for (int thisIndex = 0; thisIndex < thisSize ; ++thisIndex) {
           if (prescribedDataMgr_->is_fixed(inode,thisField,thisIndex)) {
             rhs(inode,thisIndex) = 0.;
           }
         }
       }
     }
   }
   //--------------------------------------------------------
   void ATC_Coupling::set_initial_conditions()
   {
     // set fields 
     prescribedDataMgr_->set_initial_conditions(time(),
       fields_,dot_fields_,ddot_fields_,dddot_fields_);
 
     // set (all) related data
     map<FieldName,int>::const_iterator field;
     for (field = fieldSizes_.begin(); field!=fieldSizes_.end(); field++) {
       FieldName thisField = field->first;
       int thisSize = field->second;
       DENS_MAT & rhs(rhs_[thisField].set_quantity());
       for (int inode = 0; inode < nNodes_ ; ++inode) {
         for (int thisIndex = 0; thisIndex < thisSize ; ++thisIndex) {
           rhs(inode,thisIndex) = 0.;
         }
       }
     }
   }
   //--------------------------------------------------------
   void ATC_Coupling::set_sources()
   {
     // set fields
     prescribedDataMgr_->set_sources(time(),sources_);
 
   }
   //-----------------------------------------------------------------
   // this is w_a source_a
   void ATC_Coupling::compute_sources_at_atoms(const RHS_MASK & rhsMask,
                                               const FIELDS & fields,
                                               const PhysicsModel * physicsModel,
                                               FIELD_MATS & atomicSources)
   {
     if (shpFcnMask_) {
       feEngine_->compute_source(rhsMask,
                                 fields,
                                 physicsModel,
                                 atomMaterialGroupsMask_,
                                 atomicWeightsMask_->quantity(),
                                 shpFcnMask_->quantity(),
                                 shpFcnDerivsMask_->quantity(),
                                 atomicSources);
     }
     else {
       for (FIELDS::const_iterator field = fields.begin(); 
            field != fields.end(); field++) {
         FieldName thisFieldName = field->first;
         FIELDS::const_iterator fieldItr = fields.find(thisFieldName);
         const DENS_MAT & f = (fieldItr->second).quantity();
         atomicSources[thisFieldName].reset(f.nRows(),f.nCols());
       }
     }
   }
   //-----------------------------------------------------------------
   
   void ATC_Coupling::compute_atomic_sources(const RHS_MASK & fieldMask,
                                             const FIELDS & fields,
                                             FIELDS & atomicSources)
   {
 
     for (FIELDS::const_iterator field = fields.begin();
          field != fields.end(); field++) {
       FieldName thisFieldName = field->first;
       if (is_intrinsic(thisFieldName)) {
         atomicSources[thisFieldName] = 0.;
         if (fieldMask(thisFieldName,FLUX)) {
           atomicSources[thisFieldName] = boundaryFlux_[thisFieldName];
         }
         if (fieldMask(thisFieldName,PRESCRIBED_SOURCE)) {
           atomicSources[thisFieldName] -= fluxMask_*(sources_[thisFieldName].quantity());
         } 
 
 
         // add in sources from extrinsic models
         if (fieldMask(thisFieldName,EXTRINSIC_SOURCE))
           atomicSources[thisFieldName] -= fluxMask_*(extrinsicSources_[thisFieldName].quantity());
 
       }
     }
   }
   //-----------------------------------------------------------------
   void ATC_Coupling::masked_atom_domain_rhs_tangent(
     const pair<FieldName,FieldName> row_col,
     const RHS_MASK & rhsMask,
     const FIELDS & fields, 
     SPAR_MAT & stiffness,
     const PhysicsModel * physicsModel)
   {
     if (shpFcnMask_) {
       feEngine_->compute_tangent_matrix(rhsMask, row_col,
                                         fields, physicsModel, atomMaterialGroupsMask_,
                                         atomicWeightsMask_->quantity(), shpFcnMask_->quantity(),
                                         shpFcnDerivsMask_->quantity(),stiffness);
     }
     else {
       stiffness.reset(nNodes_,nNodes_);
     }
   }
   //-----------------------------------------------------------------
   void ATC_Coupling::compute_rhs_tangent(
     const pair<FieldName,FieldName> row_col,
     const RHS_MASK & rhsMask,
     const FIELDS & fields, 
     SPAR_MAT & stiffness,
     const IntegrationDomainType integrationType,
     const PhysicsModel * physicsModel)
   {
 
-    if (integrationType  == FULL_DOMAIN_ATOMIC_QUADRATURE_SOURCE) {
-      RHS_MASK rhsMaskFE = rhsMask;
-      RHS_MASK rhsMaskMD = rhsMask; rhsMaskMD = false;
-      for (FIELDS::const_iterator field = fields.begin(); 
-           field != fields.end(); field++) {
-        FieldName thisFieldName = field->first;
-        if ( rhsMaskFE(thisFieldName,SOURCE) ) {
-          rhsMaskFE(thisFieldName,SOURCE) = false;
-          rhsMaskMD(thisFieldName,SOURCE) = true;
-        }
-      }
-      feEngine_->compute_tangent_matrix(rhsMaskFE, row_col,
-        fields , physicsModel, elementToMaterialMap_, stiffness);
-      masked_atom_domain_rhs_tangent(row_col,
-                                     rhsMaskMD,
-                                     fields,
-                                     stiffnessAtomDomain_,
-                                     physicsModel);
-      stiffness += stiffnessAtomDomain_; 
-
-    }
-    else {
       feEngine_->compute_tangent_matrix(rhsMask, row_col,
         fields , physicsModel, elementToMaterialMap_, stiffness);
-    }
     ROBIN_SURFACE_SOURCE & robinFcn = *(prescribedDataMgr_->robin_functions()); 
     feEngine_->add_robin_tangent(rhsMask, fields, time(), robinFcn, stiffness);
   }
   //-----------------------------------------------------------------
   void ATC_Coupling::compute_rhs_vector(const RHS_MASK & rhsMask,
                                         const FIELDS & fields, 
                                               FIELDS & rhs,
                                         const IntegrationDomainType domain,
                                         const PhysicsModel * physicsModel)
   {
     if (!physicsModel) physicsModel = physicsModel_;
     
 
     // compute FE contributions
     
     evaluate_rhs_integral(rhsMask,fields,rhs,domain,physicsModel);
 
     for (int n = 0; n < rhsMask.nRows(); n++) {
       FieldName thisFieldName = FieldName(n);
       if (rhsMask(thisFieldName,PRESCRIBED_SOURCE)) {
         if (is_intrinsic(thisFieldName)) {
           rhs[thisFieldName] += fluxMaskComplement_*(sources_[thisFieldName].quantity());
         } 
         else {
           rhs[thisFieldName] +=                     sources_[thisFieldName].quantity();
         }
       }
 
       // add in sources from extrinsic models
       if (rhsMask(thisFieldName,EXTRINSIC_SOURCE)) {
         if (is_intrinsic(thisFieldName)) {
           rhs[thisFieldName] += fluxMaskComplement_*(extrinsicSources_[thisFieldName].quantity());
         }
         else {
           rhs[thisFieldName] +=                     extrinsicSources_[thisFieldName].quantity();
         }
       }
       
     }
     ROBIN_SURFACE_SOURCE & robinFcn = *(prescribedDataMgr_->robin_functions()); 
     feEngine_->add_robin_fluxes(rhsMask, fields, time(), robinFcn, rhs);
   }
   //-----------------------------------------------------------------
   void ATC_Coupling::masked_atom_domain_rhs_integral(
     const Array2D<bool> & rhsMask,
     const FIELDS & fields, FIELDS & rhs,
     const PhysicsModel * physicsModel)
   {
     if (shpFcnMask_) {
       feEngine_->compute_rhs_vector(rhsMask,
                                     fields,
                                     physicsModel,
                                     atomMaterialGroupsMask_,
                                     atomicWeightsMask_->quantity(),
                                     shpFcnMask_->quantity(),
                                     shpFcnDerivsMask_->quantity(),
                                     rhs);
     }
     else {
       for (FIELDS::const_iterator field = fields.begin(); 
            field != fields.end(); field++) {
         FieldName thisFieldName = field->first;
         FIELDS::const_iterator fieldItr = fields.find(thisFieldName);
         const DENS_MAT & f = (fieldItr->second).quantity();
         (rhs[thisFieldName].set_quantity()).reset(f.nRows(),f.nCols());
       }
     }
   }
   //-----------------------------------------------------------------
   void ATC_Coupling::evaluate_rhs_integral(
     const Array2D<bool> & rhsMask,
     const FIELDS & fields, FIELDS & rhs,
     const IntegrationDomainType integrationType,
     const PhysicsModel * physicsModel)
   {
     
     if (!physicsModel) physicsModel = physicsModel_;
 
     
     if      (integrationType == FE_DOMAIN ) {
       feEngine_->compute_rhs_vector(rhsMask, 
                                     fields, 
                                     physicsModel,
                                     elementToMaterialMap_,
                                     rhs, 
                                     &(elementMask_->quantity()));
       masked_atom_domain_rhs_integral(rhsMask,
                                       fields,
                                       rhsAtomDomain_,
                                       physicsModel);
       for (FIELDS::const_iterator field = fields.begin(); 
            field != fields.end(); field++) {
         FieldName thisFieldName = field->first;
         rhs[thisFieldName] -= rhsAtomDomain_[thisFieldName].quantity();
       }
     }
     else if (integrationType == ATOM_DOMAIN) {
       
       masked_atom_domain_rhs_integral(rhsMask,
                                       fields,
                                       rhs,
                                       physicsModel);
     }
-    else if (integrationType  == FULL_DOMAIN_ATOMIC_QUADRATURE_SOURCE) {
-      RHS_MASK rhsMaskFE = rhsMask;
-      RHS_MASK rhsMaskMD = rhsMask; rhsMaskMD = false;
-      for (FIELDS::const_iterator field = fields.begin(); 
-           field != fields.end(); field++) {
-        FieldName thisFieldName = field->first;
-        if ( rhsMaskFE(thisFieldName,SOURCE) ) {
-          rhsMaskFE(thisFieldName,SOURCE) = false;
-          rhsMaskMD(thisFieldName,SOURCE) = true;
-        }
-      }
-      feEngine_->compute_rhs_vector(rhsMaskFE,
-                                    fields,
-                                    physicsModel,
-                                    elementToMaterialMap_,
-                                    rhs);
-      masked_atom_domain_rhs_integral(rhsMaskMD,
-                                      fields,
-                                      rhsAtomDomain_,
-                                      physicsModel);
-      for (FIELDS::const_iterator field = fields.begin(); 
-           field != fields.end(); field++) {
-        FieldName thisFieldName = field->first;
-
-        if ( ((rhs[thisFieldName].quantity()).size() > 0)
-         && ((rhsAtomDomain_[thisFieldName].quantity()).size() > 0) )
-          rhs[thisFieldName] += rhsAtomDomain_[thisFieldName].quantity();
-      }
-    }
     else { // domain == FULL_DOMAIN
       feEngine_->compute_rhs_vector(rhsMask,
                                     fields,
                                     physicsModel,
                                     elementToMaterialMap_,
                                     rhs);
     }
   }
   
   //--------------------------------------------------
   bool ATC_Coupling::reset_methods() const
   {
     bool resetMethods = ATC_Method::reset_methods() || atomicRegulator_->need_reset();
     for (_ctiIt_ = timeIntegrators_.begin(); _ctiIt_ != timeIntegrators_.end(); ++_ctiIt_) {
       resetMethods |= (_ctiIt_->second)->need_reset();
     }
     return resetMethods;
   }
   //--------------------------------------------------
   void ATC_Coupling::initialize()
   { 
     // initialize physics model
     if (physicsModel_) physicsModel_->initialize();
 
     ATC_Method::initialize();
     
     // initialized_ is set to true by derived class initialize()
     // STEP 6 - data initialization continued:  set initial conditions
     if (!initialized_) {
       // Apply integration masking and new ICs
       // initialize schedule derivatives
       try {
         set_initial_conditions();
       }
       catch (ATC::ATC_Error& atcError) {
         if (!useRestart_)
           throw;
       }
     }
    
     // initialize and fix computational geometry, this can be changed in the future for Eulerian calculations that fill and empty elements which is why it is outside a !initialized_ guard
     internalElement_->unfix_quantity();
     if (ghostElement_) ghostElement_->unfix_quantity();
     internalElement_->quantity();
     if (ghostElement_) ghostElement_->quantity();
     nodalGeometryType_->quantity();
     internalElement_->fix_quantity();
     if (ghostElement_) ghostElement_->fix_quantity();
     reset_flux_mask();
 
     // setup grouping of atoms by material
     reset_atom_materials();
 
     // reset time filters if needed
     if (timeFilterManager_.need_reset()) {
       if ((!initialized_) || (atomToElementMapType_ == EULERIAN)) {
         map<FieldName,int>::const_iterator field;
         for (field = fieldSizes_.begin(); field!=fieldSizes_.end(); field++) {
           FieldName thisField = field->first;
           if (is_intrinsic(thisField) && is_dynamic(thisField)) { 
             compute_mass_matrix(thisField);
             if (!useConsistentMassMatrix_(thisField) && !useFeMdMassMatrix_) {
               massMatsMd_[thisField] = massMatsMdInstantaneous_[thisField].quantity();
               massMatsAq_[thisField] = massMatsAqInstantaneous_[thisField].quantity();
               update_mass_matrix(thisField);
             }
           }
         }
       }
     }
     
     // prepare computes for first timestep
     lammpsInterface_->computes_addstep(lammpsInterface_->ntimestep()+1);
 
     // resetting precedence:
     // time integrator -> kinetostat/thermostat -> time filter
     // init_filter uses fieldRateNdFiltered which comes from the time integrator,
     // which is why the time integrator is initialized first
 
     // other initializations
     if (reset_methods()) {
       for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
         (_tiIt_->second)->initialize();
       }
       atomicRegulator_->initialize();
     }
     extrinsicModelManager_.initialize();
     if (timeFilterManager_.need_reset()) {// reset thermostat power
       init_filter();
     }
     // clears need for reset
     timeFilterManager_.initialize();
     ghostManager_.initialize();
 
     if (!initialized_) {
       // initialize sources based on initial FE temperature
       double dt = lammpsInterface_->dt();
       prescribedDataMgr_->set_sources(time()+0.5*dt,sources_);
       extrinsicModelManager_.set_sources(fields_,extrinsicSources_);
       atomicRegulator_->compute_boundary_flux(fields_);
       compute_atomic_sources(fieldMask_,fields_,atomicSources_);
 
       // read in field data if necessary
       if (useRestart_) {
         RESTART_LIST data;
         read_restart_data(restartFileName_,data);
         useRestart_ = false;
       }
 
       // set consistent initial conditions, if requested
       if (!timeFilterManager_.filter_dynamics() && consistentInitialization_) {
         
         const INT_ARRAY & nodeType(nodalGeometryType_->quantity());
 
         if (fieldSizes_.find(VELOCITY) != fieldSizes_.end()) {
           DENS_MAT & velocity(fields_[VELOCITY].set_quantity());
           DENS_MAN * nodalAtomicVelocity(interscaleManager_.dense_matrix("NodalAtomicVelocity"));
           const DENS_MAT & atomicVelocity(nodalAtomicVelocity->quantity());
           for (int i = 0; i<nNodes_; ++i) {
             
             if (nodeType(i,0)==MD_ONLY) {
               for (int j = 0; j < nsd_; j++) {
                 velocity(i,j) = atomicVelocity(i,j);
               }
             }
           }
         }
 
         if (fieldSizes_.find(TEMPERATURE) != fieldSizes_.end()) {
           DENS_MAT & temperature(fields_[TEMPERATURE].set_quantity());
           DENS_MAN * nodalAtomicTemperature(interscaleManager_.dense_matrix("NodalAtomicTemperature"));
           const DENS_MAT & atomicTemperature(nodalAtomicTemperature->quantity());
             
           for (int i = 0; i<nNodes_; ++i) {
             
             if (nodeType(i,0)==MD_ONLY) {
               temperature(i,0) = atomicTemperature(i,0);
             }
           }
         }
 
         if (fieldSizes_.find(DISPLACEMENT) != fieldSizes_.end()) {
           DENS_MAT & displacement(fields_[DISPLACEMENT].set_quantity());
           DENS_MAN * nodalAtomicDisplacement(interscaleManager_.dense_matrix("NodalAtomicDisplacement"));
           const DENS_MAT & atomicDisplacement(nodalAtomicDisplacement->quantity());
           for (int i = 0; i<nNodes_; ++i) {
             
             if (nodeType(i,0)==MD_ONLY) {
               for (int j = 0; j < nsd_; j++) {
                 displacement(i,j) = atomicDisplacement(i,j);
               }
             }
           }
         }
           
         //WIP_JAT update next two when full species time integrator is added
         if (fieldSizes_.find(MASS_DENSITY) != fieldSizes_.end()) {
           DENS_MAT & massDensity(fields_[MASS_DENSITY].set_quantity());
           const DENS_MAT & atomicMassDensity(atomicFields_[MASS_DENSITY]->quantity());
           for (int i = 0; i<nNodes_; ++i) {
             
             if (nodeType(i,0)==MD_ONLY) {
               massDensity(i,0) = atomicMassDensity(i,0);
             }
           }
         }
           
         if (fieldSizes_.find(SPECIES_CONCENTRATION) != fieldSizes_.end()) {
           DENS_MAT & speciesConcentration(fields_[SPECIES_CONCENTRATION].set_quantity());
           const DENS_MAT & atomicSpeciesConcentration(atomicFields_[SPECIES_CONCENTRATION]->quantity());
           for (int i = 0; i<nNodes_; ++i) {
             
             if (nodeType(i,0)==MD_ONLY) {
               for (int j = 0; j < atomicSpeciesConcentration.nCols(); ++j) {
                 speciesConcentration(i,j) = atomicSpeciesConcentration(i,j);
               }
             }
           }
         }
       }
       
       initialized_ = true;
     }
 
   }
   //-------------------------------------------------------------------
   void ATC_Coupling::construct_time_integration_data()
   {
     if (!initialized_) {
 
       map<FieldName,int>::const_iterator field;
       for (field = fieldSizes_.begin(); field!=fieldSizes_.end(); field++) {
         FieldName thisField = field->first;
         int thisSize = field->second;
         
         // Allocate fields, initialize to default values, set up initial schedule
         
         fields_[thisField].reset(nNodes_,thisSize);
         dot_fields_[thisField].reset(nNodes_,thisSize);
         ddot_fields_[thisField].reset(nNodes_,thisSize);
         dddot_fields_[thisField].reset(nNodes_,thisSize);
         
         // Allocate restricted fields
         if (is_intrinsic(thisField)) {
           nodalAtomicFields_[thisField].reset(nNodes_,thisSize);
           nodalAtomicFieldsRoc_[thisField].reset(nNodes_,thisSize);
         }
 
         // Dimension finite element rhs matrix
         rhs_[thisField].reset(nNodes_,thisSize);
         rhsAtomDomain_[thisField].reset(nNodes_,thisSize);
         
         sources_[thisField].reset(nNodes_,thisSize);
         extrinsicSources_[thisField].reset(nNodes_,thisSize);
         boundaryFlux_[thisField].reset(nNodes_,thisSize);
         
         if (is_intrinsic(thisField) && is_dynamic(thisField)) {
           massMats_[thisField].reset(nNodes_,nNodes_); // PARALLELIZE
           massMatsFE_[thisField].reset(nNodes_,nNodes_);
           massMatsAq_[thisField].reset(nNodes_,nNodes_);
           massMatsMd_[thisField].reset(nNodes_,nNodes_);
           massMatsMdInstantaneous_[thisField].reset(nNodes_,nNodes_);
           massMatsAqInstantaneous_[thisField].reset(nNodes_,nNodes_);
           massMatsInv_[thisField].reset(nNodes_,nNodes_);  // PARALLELIZE
           massMatsMdInv_[thisField].reset(nNodes_,nNodes_); // PARALLELIZE
         }
         else {
           // no MD mass matrices needed, regular matrices computed in extrinsic model
           if (useConsistentMassMatrix_(thisField)) {
             // compute FE mass matrix in full domain
             
             consistentMassMats_[thisField].reset(nNodes_,nNodes_); // PARALLELIZE
             consistentMassMatsInv_[thisField].reset(nNodes_,nNodes_); // PARALLELIZE
           }
           else {
             massMats_[thisField].reset(nNodes_,nNodes_); // PARALLELIZE
             massMatsInv_[thisField].reset(nNodes_,nNodes_); // PARALLELIZE
           }
         }
       }
     }
   }
   //--------------------------------------------------------
   //  create_full_element_mask
   //    constructs element mask which only masks out 
   //    null elements
   //--------------------------------------------------------
   MatrixDependencyManager<DenseMatrix, bool> * ATC_Coupling::create_full_element_mask()
   {
     MatrixDependencyManager<DenseMatrix, bool> * elementMaskMan = new MatrixDependencyManager<DenseMatrix, bool>(feEngine_->num_elements(),1);
     DenseMatrix<bool> & elementMask(elementMaskMan->set_quantity());
     elementMask = true;
       
     const set<int> & nullElements = feEngine_->null_elements();
     set<int>::const_iterator iset;
     for (iset = nullElements.begin(); iset != nullElements.end(); iset++) {
       int ielem = *iset;
       elementMask(ielem,0) = false;
     }
 
     return elementMaskMan;
   }
   //--------------------------------------------------------
   //  create_element_set_mask
   //    constructs element mask based on an element set,
   //    uses ints for MPI communication later
   //--------------------------------------------------------
   MatrixDependencyManager<DenseMatrix, int> * ATC_Coupling::create_element_set_mask(const string & elementSetName)
   {
     MatrixDependencyManager<DenseMatrix, int> * elementMaskMan = new MatrixDependencyManager<DenseMatrix, int>(feEngine_->num_elements(),1);
     DenseMatrix<int> & elementMask(elementMaskMan->set_quantity());
     elementMask = false;
 
     const set<int> & elementSet((feEngine_->fe_mesh())->elementset(elementSetName));
     set<int>::const_iterator iset;
     for (iset = elementSet.begin(); iset != elementSet.end(); ++iset) {
       int ielem = *iset;
       elementMask(ielem,0) = true;
     }
       
     const set<int> & nullElements = feEngine_->null_elements();
     for (iset = nullElements.begin(); iset != nullElements.end(); iset++) {
       int ielem = *iset;
       elementMask(ielem,0) = false;
     }
 
     return elementMaskMan;
   }
   //--------------------------------------------------------
   //  set_computational_geometry
   //    constructs needed transfer operators which define
   //    hybrid atom/FE computational geometry
   //--------------------------------------------------------
   void ATC_Coupling::set_computational_geometry()
   {
     ATC_Method::set_computational_geometry();
 
     // does element contain internal atoms
     if (internalElementSet_.size()) {
       // set up elements and maps based on prescribed element sets
       internalElement_ = create_element_set_mask(internalElementSet_);
     }
     else {
       internalElement_ = new AtomTypeElement(this,atomElement_);
     }
     interscaleManager_.add_dense_matrix_int(internalElement_,
                                             "ElementHasInternal");
 
     if (groupbitGhost_) {
       atomGhostElement_ = new AtomToElementMap(this,
                                                atomGhostCoarseGrainingPositions_,
                                                GHOST);
       interscaleManager_.add_per_atom_int_quantity(atomGhostElement_,
                                                    "AtomGhostElement");
       
       // does element contain ghost atoms
       ghostElement_ = new AtomTypeElement(this,atomGhostElement_);
       interscaleManager_.add_dense_matrix_int(ghostElement_,
                                               "ElementHasGhost");
     }
     
     // element masking for approximate right-hand side FE atomic quadrature
     if (atomQuadForInternal_) {
       elementMask_ = create_full_element_mask();
     }
     else {
       if (internalElementSet_.size()) {
         // when geometry is based on elements, there are no mixed elements
         elementMask_ = new MatrixDependencyManager<DenseMatrix, bool>;
         (elementMask_->set_quantity()).reset(feEngine_->num_elements(),1,false);
       }
       else {
         elementMask_ = new ElementMask(this);
       }
       internalToMask_ = new AtomToElementset(this,elementMask_);
       interscaleManager_.add_per_atom_int_quantity(internalToMask_,
                                                    "InternalToMaskMap");
     }
     interscaleManager_.add_dense_matrix_bool(elementMask_,
                                              "ElementMask");
 
     if (useFeMdMassMatrix_) {
       if (atomQuadForInternal_) {
         elementMaskMass_ = elementMask_;
       }
       else {
         elementMaskMass_ = create_full_element_mask();
         interscaleManager_.add_dense_matrix_bool(elementMaskMass_,
                                                  "NonNullElementMask");
       }
 
       elementMaskMassMd_ = new AtomElementMask(this);
       interscaleManager_.add_dense_matrix_bool(elementMaskMassMd_,
                                                "InternalElementMask");
     }
 
     // assign element and node types for computational geometry
     if (internalElementSet_.size()) {
       nodalGeometryType_ = new NodalGeometryTypeElementSet(this);
     }
     else {
       nodalGeometryType_ = new NodalGeometryType(this);
     }
     interscaleManager_.add_dense_matrix_int(nodalGeometryType_,
                                             "NodalGeometryType");
   }
   //--------------------------------------------------------
   //  construct_interpolant
   //    constructs: interpolatn, accumulant, weights, and spatial derivatives
   //--------------------------------------------------------
   void ATC_Coupling::construct_interpolant()
   {
     // finite element shape functions for interpolants
     PerAtomShapeFunction * atomShapeFunctions = new PerAtomShapeFunction(this);
     interscaleManager_.add_per_atom_sparse_matrix(atomShapeFunctions,"Interpolant");
     shpFcn_ = atomShapeFunctions;
 
     // use shape functions for accumulants if no kernel function is provided
     if (!kernelFunction_) {
       accumulant_ = shpFcn_;
     }
     else {
       if (kernelOnTheFly_) throw ATC_Error("ATC_Coupling::construct_transfers - on the fly kernel evaluations not supported");
       PerAtomKernelFunction * atomKernelFunctions = new PerAtomKernelFunction(this);
       interscaleManager_.add_per_atom_sparse_matrix(atomKernelFunctions,
                                                     "Accumulant");
       accumulant_ = atomKernelFunctions;
       accumulantWeights_ = new AccumulantWeights(accumulant_);
       mdMassNormalization_ = false;
     }
     
     this->create_atom_volume();
 
     // masked atom weights
     if (atomQuadForInternal_) {
       atomicWeightsMask_ = atomVolume_;
     }
     else {
       atomicWeightsMask_ = new MappedDiagonalMatrix(this,
                                                     atomVolume_,
                                                     internalToMask_);
       interscaleManager_.add_diagonal_matrix(atomicWeightsMask_,
                                              "AtomWeightsMask");
     }
     // nodal volumes for mass matrix, relies on atomVolumes constructed in base class construct_transfers
     nodalAtomicVolume_ = new AdmtfShapeFunctionRestriction(this,atomVolume_,shpFcn_);
     interscaleManager_.add_dense_matrix(nodalAtomicVolume_,"NodalAtomicVolume");
 
     // shape function derivatives, masked shape function and derivatives if needed for FE quadrature in atomic domain
     if (atomQuadForInternal_) {
       shpFcnDerivs_ = new PerAtomShapeFunctionGradient(this);
       interscaleManager_.add_vector_sparse_matrix(shpFcnDerivs_,
                                                   "InterpolantGradient");
 
       shpFcnMask_ = shpFcn_;
       shpFcnDerivsMask_ = shpFcnDerivs_;
     }
     else {
       bool hasMaskedElt = false;
       const DenseMatrix<bool> & elementMask(elementMask_->quantity());
       for (int i = 0; i < elementMask.size(); ++i) {
         if (elementMask(i,0)) {
           hasMaskedElt = true;
           break;
         }
       }
       if (hasMaskedElt) {
         shpFcnDerivs_ = new PerAtomShapeFunctionGradient(this);
         interscaleManager_.add_vector_sparse_matrix(shpFcnDerivs_,
                                                     "InterpolantGradient");
 
         shpFcnMask_ = new RowMappedSparseMatrix(this,
                                                 shpFcn_,
                                                 internalToMask_);
         interscaleManager_.add_sparse_matrix(shpFcnMask_,
                                              "ShapeFunctionMask");
         shpFcnDerivsMask_ = new RowMappedSparseMatrixVector(this,
                                                             shpFcnDerivs_,
                                                             internalToMask_);
         interscaleManager_.add_vector_sparse_matrix(shpFcnDerivsMask_,"ShapeFunctionGradientMask");
       }
     }
   }
   //--------------------------------------------------------
   //  construct_molecule_transfers
   //--------------------------------------------------------
   void ATC_Coupling::construct_molecule_transfers()
   {
     
     map<string,pair<MolSize,int> >::const_iterator molecule;
     PerAtomQuantity<double> * atomProcGhostCoarseGrainingPositions = interscaleManager_.per_atom_quantity("AtomicProcGhostCoarseGrainingPositions");
     FundamentalAtomQuantity * mass = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_MASS,
                                                                                   PROC_GHOST);
     for (molecule = moleculeIds_.begin(); molecule != moleculeIds_.end(); molecule++) {
       const string moleculeName = molecule->first;
       int groupbit = (molecule->second).second;
       SmallMoleculeSet * smallMoleculeSet = new SmallMoleculeSet(this,groupbit);
       smallMoleculeSet->initialize();
       interscaleManager_.add_small_molecule_set(smallMoleculeSet,moleculeName);
       SmallMoleculeCentroid * moleculeCentroid = 
         new SmallMoleculeCentroid(this,mass,smallMoleculeSet,atomProcGhostCoarseGrainingPositions);
       interscaleManager_.add_dense_matrix(moleculeCentroid,"MoleculeCentroid"+moleculeName);
 
       // shape function at molecular coordinates
       PointToElementMap * elementMapMol = 
         new PointToElementMap(this,moleculeCentroid);
       interscaleManager_.add_dense_matrix_int(elementMapMol,
                                               "ElementMap"+moleculeName);
       InterpolantSmallMolecule * shpFcnMol = new InterpolantSmallMolecule(this,
         elementMapMol, moleculeCentroid, smallMoleculeSet);
       interscaleManager_.add_sparse_matrix(shpFcnMol,
                                            "ShapeFunction"+moleculeName);
     }
   }
   //--------------------------------------------------------
   //  construct_transfers
   //    constructs needed transfer operators
   //--------------------------------------------------------
   void ATC_Coupling::construct_transfers()
   {
     ATC_Method::construct_transfers();
 
     if (!useFeMdMassMatrix_) {
       // transfer for MD mass matrices based on requested intrinsic fields
       if (fieldSizes_.find(TEMPERATURE) != fieldSizes_.end()) {
         // classical thermodynamic heat capacity of the atoms
         HeatCapacity * heatCapacity = new HeatCapacity(this);
         interscaleManager_.add_per_atom_quantity(heatCapacity,
                                                  "AtomicHeatCapacity");
 
         // atomic thermal mass matrix
         nodalAtomicHeatCapacity_ = new AtfShapeFunctionRestriction(this,
                                                                    heatCapacity,
                                                                    shpFcn_);
         interscaleManager_.add_dense_matrix(nodalAtomicHeatCapacity_,
                                             "NodalAtomicHeatCapacity");
       }
       if ((fieldSizes_.find(VELOCITY) != fieldSizes_.end()) || (fieldSizes_.find(DISPLACEMENT) != fieldSizes_.end())) {
         // atomic momentum mass matrix
         FundamentalAtomQuantity * atomicMass = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_MASS);
         nodalAtomicMass_ = new AtfShapeFunctionRestriction(this,
                                                            atomicMass,
                                                            shpFcn_);
         interscaleManager_.add_dense_matrix(nodalAtomicMass_,
                                             "AtomicMomentumMassMat");
       }
       if (fieldSizes_.find(MASS_DENSITY) != fieldSizes_.end()) {
         // atomic dimensionless mass matrix
         ConstantQuantity<double> * atomicOnes = new ConstantQuantity<double>(this,1);
         interscaleManager_.add_per_atom_quantity(atomicOnes,"AtomicOnes");
         nodalAtomicCount_ = new AtfShapeFunctionRestriction(this,
                                                             atomicOnes,
                                                             shpFcn_);
         interscaleManager_.add_dense_matrix(nodalAtomicCount_,
                                             "AtomicDimensionlessMassMat");
       }
     }
 
     extrinsicModelManager_.construct_transfers();
   }
   //--------------------------------------------------
   void ATC_Coupling::delete_mass_mat_time_filter(FieldName thisField)
   {
   }
   //--------------------------------------------------
   void ATC_Coupling::set_mass_mat_time_filter(FieldName thisField,TimeFilterManager::FilterIntegrationType filterIntegrationType)
   {
     massMatTimeFilters_[thisField] = timeFilterManager_.construct(filterIntegrationType);
   }
   //--------------------------------------------------------------
   /** method to trigger construction of mesh data after mesh construction */
   //--------------------------------------------------------------
   void ATC_Coupling::initialize_mesh_data(void)
   {
     int nelts = feEngine_->fe_mesh()->num_elements();
     elementToMaterialMap_.reset(nelts);
     elementToMaterialMap_ = 0;
 
     construct_prescribed_data_manager();
     meshDataInitialized_ = true;
   }
   //--------------------------------------------------------
   
   void ATC_Coupling::reset_flux_mask(void)
   {
     int i;
     // this is exact only for uniform meshes and certain types of atomic weights
     // \int_{\Omega_MD} N_I dV = \sum_\alpha N_I\alpha V_\alpha
     fluxMask_.reset((invNodeVolumes_.quantity()) 
               * (nodalAtomicVolume_->quantity()));
 
     DIAG_MAT id(fluxMask_.nRows(),fluxMask_.nCols()); 
     id = 1.0;
     fluxMaskComplement_ = id + -1.0*fluxMask_;
 
     // set flux masks for nodes we can tell by geometry
     const INT_ARRAY & nodeType(nodalGeometryType_->quantity());
     for (i = 0; i < nNodes_; ++i) {
       if (nodeType(i,0)==MD_ONLY) {
         fluxMask_(i,i) = 1.;
         fluxMaskComplement_(i,i) = 0.;
       }
       else if (nodeType(i,0)==FE_ONLY) {
         fluxMask_(i,i) = 0.;
         fluxMaskComplement_(i,i) = 1.;
       }
     }
   }
 
   //--------------------------------------------------------
   void ATC_Coupling::compute_mass_matrix(FieldName thisField, PhysicsModel * physicsModel)
   {
 
     if (!physicsModel) physicsModel = physicsModel_;
     if (useConsistentMassMatrix_(thisField)) {
       // compute FE mass matrix in full domain
       
       Array<FieldName> massMask(1);
       massMask(0) = thisField; 
       
       feEngine_->compute_mass_matrix(massMask,fields_,physicsModel,
                                      elementToMaterialMap_,consistentMassMats_,
                                      &(elementMask_->quantity()));
       // brute force computation of inverse
       consistentMassMatsInv_[thisField] = inv((consistentMassMats_[thisField].quantity()).dense_copy());
     }
     else { // lumped mass matrix
       // compute FE mass matrix in full domain
       Array<FieldName> massMask(1);
       massMask(0) = thisField;
 
       if (useFeMdMassMatrix_) {
         feEngine_->compute_lumped_mass_matrix(massMask,fields_,physicsModel,
                                               elementToMaterialMap_,massMats_,
                                               &(elementMaskMass_->quantity()));
         const DIAG_MAT & myMassMat(massMats_[thisField].quantity());
         DIAG_MAT & myMassMatInv(massMatsInv_[thisField].set_quantity());
         DIAG_MAT & myMassMatMdInv(massMatsMdInv_[thisField].set_quantity());
 
         feEngine_->compute_lumped_mass_matrix(massMask,fields_,physicsModel,
                                               elementToMaterialMap_,massMatsMd_,
                                               &(elementMaskMassMd_->quantity()));
         const DIAG_MAT & myMassMatMd(massMatsMd_[thisField].quantity());
         // compute inverse mass matrices since we're using lumped masses
         for (int iNode = 0; iNode < nNodes_; iNode++) {
           
           if (fabs(myMassMat(iNode,iNode))>0)
             myMassMatInv(iNode,iNode) = 1./myMassMat(iNode,iNode);
           else
             myMassMatInv(iNode,iNode) = 0.;
 
           if (fabs(myMassMatMd(iNode,iNode))>0)
             myMassMatMdInv(iNode,iNode) = 1./myMassMatMd(iNode,iNode);
           else
             myMassMatMdInv(iNode,iNode) = 0.;
         }
       }
       else {
         feEngine_->compute_lumped_mass_matrix(massMask,fields_,physicsModel,
                                               elementToMaterialMap_,massMatsFE_,
                                               &(elementMask_->quantity()));
         // fully remove contributions from internal nodes
         
         DIAG_MAT & myMassMatFE(massMatsFE_[thisField].set_quantity());
         if (!atomQuadForInternal_) {
           const INT_ARRAY & nodeType(nodalGeometryType_->quantity());
           for (int iNode = 0; iNode < nNodes_; iNode++)
             if (nodeType(iNode,0)==MD_ONLY) {
               myMassMatFE(iNode,iNode) = 0.;
             }
         }
         
         // atomic quadrature for FE mass matrix in atomic domain
         if (shpFcnMask_) {
           feEngine_->compute_lumped_mass_matrix(massMask,fields_,physicsModel,atomMaterialGroupsMask_,
                                                 atomicWeightsMask_->quantity(),shpFcnMask_->quantity(),
                                                 massMatsAqInstantaneous_);
         }
         else {
           (massMatsAqInstantaneous_[thisField].set_quantity()).reset(nNodes_,nNodes_);
         }
         
         // set up mass MD matrices
         compute_md_mass_matrix(thisField,massMatsMdInstantaneous_[thisField].set_quantity());
       }
     }
   }
   //--------------------------------------------------------
   void ATC_Coupling::update_mass_matrix(FieldName thisField)
   {
     DIAG_MAT & myMassMat(massMats_[thisField].set_quantity());
     DIAG_MAT & myMassMatInv(massMatsInv_[thisField].set_quantity());
     DIAG_MAT & myMassMatMDInv(massMatsMdInv_[thisField].set_quantity());
     const DIAG_MAT & myMassMatMD(massMatsMd_[thisField].quantity());
 
     myMassMat = massMatsFE_[thisField].quantity();
     // remove contributions from overlap by approximate quadrature
     myMassMat -= massMatsAq_[thisField].quantity();
     // add contributions from atomic region
     myMassMat += myMassMatMD;
 
     // compute inverse mass matrices since we're using lumped masses
     for (int iNode = 0; iNode < nNodes_; iNode++) {
       
       if (fabs(myMassMatMD(iNode,iNode))>0) { 
         myMassMatMDInv(iNode,iNode) = 1./myMassMatMD(iNode,iNode);
       }
       else
         myMassMatMDInv(iNode,iNode) = 0.;
       
       if (fabs(myMassMat(iNode,iNode))>0) {
         myMassMatInv(iNode,iNode) = 1./myMassMat(iNode,iNode);
       }
       else
         myMassMatInv(iNode,iNode) = 0.;
     }
   }
 
   //---------------------------------------------------------
   //  compute_md_mass_matrix
   //    compute the mass matrix arising from only atomistic
   //    quadrature and contributions as a summation
   //---------------------------------------------------------
   void ATC_Coupling::compute_md_mass_matrix(FieldName thisField,
                                             DIAG_MAT & massMat)
   {
     
     if (thisField == TEMPERATURE) {
       massMat.shallowreset(nodalAtomicHeatCapacity_->quantity());
     }
     
     else if (thisField == DISPLACEMENT || thisField == VELOCITY) {
       massMat.shallowreset(nodalAtomicMass_->quantity());
     }
     else if (thisField == MASS_DENSITY || thisField == SPECIES_CONCENTRATION) {
       massMat.shallowreset(nodalAtomicVolume_->quantity());
     }
   }   
 
   //--------------------------------------------------
   // write_restart_file
   //   bundle matrices that need to be saved and call
   //   fe_engine to write the file
   //--------------------------------------------------
   void ATC_Coupling::write_restart_data(string fileName, RESTART_LIST & data)
   {
     atomicRegulator_->pack_fields(data);
     ATC_Method::write_restart_data(fileName,data);
   }
   
   //--------------------------------------------------
   // read_restart_file
   //   bundle matrices that need to be saved and call
   //   fe_engine to write the file
   //--------------------------------------------------
   void ATC_Coupling::read_restart_data(string fileName, RESTART_LIST & data)
   {
     atomicRegulator_->pack_fields(data);
     ATC_Method::read_restart_data(fileName,data);
   }
 
   //--------------------------------------------------
   void ATC_Coupling::reset_nlocal()
   {
     ATC_Method::reset_nlocal();
     atomicRegulator_->reset_nlocal();
   }
 
   //--------------------------------------------------------
   void ATC_Coupling::reset_atom_materials()
   {
     int nMaterials = physicsModel_->nMaterials();
     atomMaterialGroups_.reset(nMaterials);
     atomMaterialGroupsMask_.reset(nMaterials);
 
     for (int i = 0; i < nMaterials; i++) {
       atomMaterialGroups_(i).clear();
       atomMaterialGroupsMask_(i).clear();
     }
 
     const INT_ARRAY & atomToElementMap(atomElement_->quantity());
     for (int i = 0; i < nLocal_; i++) {
       atomMaterialGroups_(elementToMaterialMap_(atomToElementMap(i,0))).insert(i);
     }
     if (atomQuadForInternal_) {
       for (int i = 0; i < nLocal_; i++) {
         atomMaterialGroupsMask_(elementToMaterialMap_(atomToElementMap(i,0))).insert(i);
       }
     }
     else {
       const INT_ARRAY & map(internalToMask_->quantity());
       for (int i = 0; i < nLocal_; i++) {
         int idx = map(i,0);
         if (idx > -1) {
           atomMaterialGroupsMask_(elementToMaterialMap_(atomToElementMap(i,0))).insert(idx);
         }
       }
     }
 
     atomicRegulator_->reset_atom_materials(elementToMaterialMap_,
                                            atomElement_);
   }
 
   //--------------------------------------------------------
   //  pre_init_integrate
   //    time integration before the lammps atomic
   //    integration of the Verlet step 1
   //--------------------------------------------------------
   void ATC_Coupling::pre_init_integrate()
   {
     ATC_Method::pre_init_integrate();
     double dt = lammpsInterface_->dt();
 
     // Perform any initialization, no actual integration
     for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
       (_tiIt_->second)->pre_initial_integrate1(dt);
     }
 
     // Apply controllers to atom velocities, if needed
     atomicRegulator_->apply_pre_predictor(dt,lammpsInterface_->ntimestep());
 
     // predict nodal fields and time derivatives
     for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
       (_tiIt_->second)->pre_initial_integrate2(dt);
     }
     extrinsicModelManager_.pre_init_integrate();
   }
 
   //--------------------------------------------------------
   //  init_integrate
   //    time integration of lammps atomic quantities
   //--------------------------------------------------------
   void ATC_Coupling::init_integrate()
   {
     atomTimeIntegrator_->init_integrate_velocity(dt());
     ghostManager_.init_integrate_velocity(dt());
     // account for other fixes doing time integration
     interscaleManager_.fundamental_force_reset(LammpsInterface::ATOM_VELOCITY);
 
     // apply constraints to velocity
     atomicRegulator_->apply_mid_predictor(dt(),lammpsInterface_->ntimestep());
 
     atomTimeIntegrator_->init_integrate_position(dt());
     ghostManager_.init_integrate_position(dt());
     // account for other fixes doing time integration
     interscaleManager_.fundamental_force_reset(LammpsInterface::ATOM_POSITION);
   }
 
   ///--------------------------------------------------------
   //  post_init_integrate
   //    time integration after the lammps atomic updates of
   //    Verlet step 1
   //--------------------------------------------------------
   void ATC_Coupling::post_init_integrate()
   {
     double dt = lammpsInterface_->dt();
   
     // Compute nodal velocity at n+1
     for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
       (_tiIt_->second)->post_initial_integrate1(dt);
     }
 
     // Update kinetostat quantities if displacement is being regulated
     atomicRegulator_->apply_post_predictor(dt,lammpsInterface_->ntimestep());
 
     // Update extrisic model
     extrinsicModelManager_.post_init_integrate();
 
     // fixed values, non-group bcs handled through FE
     set_fixed_nodes();
       
     update_time(0.5);
 
     // ghost update, if needed
     ATC_Method::post_init_integrate();
 
     // Apply time filtering to mass matrices, if needed
     if ((atomToElementMapType_ == EULERIAN) && timeFilterManager_.filter_dynamics() && !useFeMdMassMatrix_) {
       map<FieldName,int>::const_iterator field;
       for (field = fieldSizes_.begin(); field!=fieldSizes_.end(); field++) {
         FieldName thisField = field->first;
         if (!useConsistentMassMatrix_(thisField) && is_intrinsic(thisField)) {
           massMatTimeFilters_[thisField]->apply_pre_step1(massMatsAq_[thisField].set_quantity(),
                                                           massMatsAqInstantaneous_[thisField].quantity(),dt);
           massMatTimeFilters_[thisField]->apply_pre_step1(massMatsMd_[thisField].set_quantity(),
                                                           massMatsMdInstantaneous_[thisField].quantity(),dt);
         }
       }
     }
   }
 
   
   //--------------------------------------------------------
   void ATC_Coupling::pre_neighbor()
   {
     ATC_Method::pre_neighbor();
     reset_atom_materials();
   }
 
   //--------------------------------------------------------
   void ATC_Coupling::pre_exchange()
   {
     ATC_Method::pre_exchange();
   }
 
   //--------------------------------------------------------
   //  pre_force
   //    prior to calculation of forces
   //--------------------------------------------------------
   void ATC_Coupling::pre_force()
   {
     ATC_Method::pre_force();
     atomicRegulator_->pre_force(); 
   }
 
   //--------------------------------------------------------
   void ATC_Coupling::post_force()
   {
     ATC_Method::post_force();
 
     if ( (atomToElementMapType_ == EULERIAN) && (step() % atomToElementMapFrequency_ == 0) ) {
       reset_atom_materials();
 
       if (!useFeMdMassMatrix_) {
         map<FieldName,int>::const_iterator field;
         for (field = fieldSizes_.begin(); field!=fieldSizes_.end(); field++) {
           FieldName thisField = field->first;
           if (is_intrinsic(thisField) && is_dynamic(thisField)) {
             compute_mass_matrix(thisField);  
           }
         }
       }
     }
 
     if (atomToElementMapType_ == EULERIAN && !useFeMdMassMatrix_) {
       if (timeFilterManager_.filter_dynamics() || (step() % atomToElementMapFrequency_ == 0)) {
         double dt = lammpsInterface_->dt();
         map<FieldName,int>::const_iterator field;
         for (field = fieldSizes_.begin(); field!=fieldSizes_.end(); field++) {
           FieldName thisField = field->first;
           if (is_intrinsic(thisField) && is_dynamic(thisField)) {
             massMatTimeFilters_[thisField]->apply_post_step1(massMatsAq_[thisField].set_quantity(),
                                                              massMatsAqInstantaneous_[thisField].quantity(),dt);
             massMatTimeFilters_[thisField]->apply_post_step1(massMatsMd_[thisField].set_quantity(),
                                                              massMatsMdInstantaneous_[thisField].quantity(),dt);
             update_mass_matrix(thisField); 
           }
         }
       }
     }
 
     // apply extrinsic model
     extrinsicModelManager_.post_force();
   }
 
   //--------------------------------------------------------
   //  post_final_integrate
   //    integration after the second stage lammps atomic 
   //    update of Verlet step 2
   //--------------------------------------------------------
   void ATC_Coupling::post_final_integrate()
   {
     double dt = lammpsInterface_->dt();
 
     // update of atomic contributions for fractional step methods
     for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
       (_tiIt_->second)->pre_final_integrate1(dt);
     }
 
     // Set sources
     prescribedDataMgr_->set_sources(time()+0.5*dt,sources_);
     extrinsicModelManager_.pre_final_integrate();
     
     bool needsSources = false;
     for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
       if ((_tiIt_->second)->has_final_predictor()) {
         needsSources = true;
         break;
       }
     }
     if (needsSources) {
       extrinsicModelManager_.set_sources(fields_,extrinsicSources_);
       atomicRegulator_->compute_boundary_flux(fields_);
       compute_atomic_sources(intrinsicMask_,fields_,atomicSources_);
     }
     atomicRegulator_->apply_pre_corrector(dt,lammpsInterface_->ntimestep());
 
     // Compute atom-integrated rhs
     // parallel communication happens within FE_Engine
     compute_rhs_vector(intrinsicMask_,fields_,rhs_,FE_DOMAIN);
     for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
       (_tiIt_->second)->add_to_rhs();
     }
     atomicRegulator_->add_to_rhs(rhs_);
     
     // Compute and add atomic contributions to FE equations
     for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
       (_tiIt_->second)->post_final_integrate1(dt);
     }
     
    // fix nodes, non-group bcs applied through FE
     set_fixed_nodes();
         
     // corrector step extrinsic model
     extrinsicModelManager_.post_final_integrate();
 
     // set state-based sources
     needsSources = false;
     for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
       if ((_tiIt_->second)->has_final_corrector()) {
         needsSources = true;
         break;
       }
     }
     if (needsSources) {
       extrinsicModelManager_.set_sources(fields_,extrinsicSources_);
       atomicRegulator_->compute_boundary_flux(fields_);
       compute_atomic_sources(intrinsicMask_,fields_,atomicSources_);
     }
 
     // Finish update of FE velocity
     for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
       (_tiIt_->second)->post_final_integrate2(dt);
     }
     
     // apply corrector phase of thermostat
     atomicRegulator_->apply_post_corrector(dt,lammpsInterface_->ntimestep());
 
     // final phase of time integration
     for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
       (_tiIt_->second)->post_final_integrate3(dt);
     }
 
     // Fix nodes, non-group bcs applied through FE
     set_fixed_nodes();
     
     update_time(0.5);
     
     output();
     lammpsInterface_->computes_addstep(lammpsInterface_->ntimestep()+1); // adds next step to computes
     //ATC_Method::post_final_integrate();
   }
 
   //=================================================================
   //
   //=================================================================
   void ATC_Coupling::finish()
   {
     ATC_Method::finish();
     // Time integrator
     for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
       (_tiIt_->second)->finish();
     }
     atomicRegulator_->finish();
   }
 
 
 
   //=================================================================
   //
   //=================================================================
   void ATC_Coupling::compute_boundary_flux(const Array2D<bool> & rhsMask,
                                            const FIELDS & fields, 
                                            FIELDS & rhs,
                                            const Array< set <int> > atomMaterialGroups,
                                            const VectorDependencyManager<SPAR_MAT * > * shpFcnDerivs,
                                            const SPAR_MAN * shpFcn,
                                            const DIAG_MAN * atomicWeights,
                                            
                                            const MatrixDependencyManager<DenseMatrix, bool> * elementMask,
                                            const RegulatedNodes * nodeSet)
   {
     if (bndyIntType_ == FE_QUADRATURE) {
       feEngine_->compute_boundary_flux(rhsMask,
                                        fields,
                                        physicsModel_,
                                        elementToMaterialMap_,
                                        (* bndyFaceSet_),
                                        rhs);
     }
     else if (bndyIntType_ == FE_INTERPOLATION) {
       if (elementMask) {
         feEngine_->compute_boundary_flux(rhsMask,
                                          fields,
                                          physicsModel_,
                                          elementToMaterialMap_,
                                          atomMaterialGroups,
                                          atomicWeights->quantity(),
                                          shpFcn->quantity(),
                                          shpFcnDerivs->quantity(),
                                          fluxMask_,
                                          rhs,
                                          &elementMask->quantity(),
                                          &nodeSet->quantity());
       }
       else {
         feEngine_->compute_boundary_flux(rhsMask,
                                          fields,
                                          physicsModel_,
                                          elementToMaterialMap_,
                                          atomMaterialGroups_,
                                          atomVolume_->quantity(),
                                          shpFcn_->quantity(),
                                          shpFcnDerivs_->quantity(),
                                          fluxMask_,
                                          rhs);
       }
     }
     else if (bndyIntType_ == NO_QUADRATURE) {
       FIELDS::const_iterator field;
       for (field = fields.begin(); field != fields.end(); field++) {
         FieldName thisFieldName = field->first;
 
 if (thisFieldName >= rhsMask.nRows()) break;
         if (rhsMask(thisFieldName,FLUX)) {
           int nrows = (field->second).nRows();
           int ncols = (field->second).nCols();
           rhs[thisFieldName].reset(nrows,ncols);
         }
       }
     }
   }
 
   //-----------------------------------------------------------------
   void ATC_Coupling::compute_flux(const Array2D<bool> & rhsMask,
                                   const FIELDS & fields, 
                                   GRAD_FIELD_MATS & flux,
                                   const PhysicsModel * physicsModel) 
   {
     if (! physicsModel) { physicsModel = physicsModel_; }
     feEngine_->compute_flux(rhsMask,
                             fields,
                             physicsModel,
                             elementToMaterialMap_,
                             flux);
   }
 
 
   //--------------------------------------------------------
 
   void ATC_Coupling::nodal_projection(const FieldName & fieldName,
                                       const PhysicsModel * physicsModel,
                                       FIELD & field)
   {
     FIELDS rhs;
     rhs[fieldName].reset(nNodes_,field.nCols());
     Array2D <bool> rhsMask(NUM_FIELDS,NUM_FLUX);
     rhsMask = false;
     rhsMask(fieldName,SOURCE) = true;
     compute_rhs_vector(rhsMask, fields_, rhs, sourceIntegration_, physicsModel);
     const DENS_MAT & B(rhs[fieldName].quantity());
 
     field = (invNodeVolumes_.quantity())*B;
   }
 
   // parse_boundary_integration
   //   parses the boundary integration to determine
   //   the type of boundary integration being used
   //--------------------------------------------------
   
   
   BoundaryIntegrationType ATC_Coupling::parse_boundary_integration(int narg,
                                                                    char **arg,
                                                                    const set< pair<int,int> > * boundaryFaceSet)
   {
     
     int argIndex = 0;
     BoundaryIntegrationType myBoundaryIntegrationType = FE_INTERPOLATION;// default
       if (narg > 0) {
         if(strcmp(arg[argIndex],"faceset")==0) {
           argIndex++;
           myBoundaryIntegrationType = FE_QUADRATURE;
           string name(arg[argIndex]);
           boundaryFaceSet = & ( (feEngine_->fe_mesh())->faceset(name));
           set_boundary_face_set(boundaryFaceSet);
         } 
         else if (strcmp(arg[argIndex],"interpolate")==0) {
           myBoundaryIntegrationType = FE_INTERPOLATION;
         }
         else if (strcmp(arg[argIndex],"no_boundary")==0) { 
           myBoundaryIntegrationType = NO_QUADRATURE; 
         }
         else {
           throw ATC_Error("Bad boundary integration type");
         }
       }
     set_boundary_integration_type(myBoundaryIntegrationType);
     return myBoundaryIntegrationType;
   }
 
 }; // namespace ATC
diff --git a/lib/atc/ATC_Coupling.h b/lib/atc/ATC_Coupling.h
index d9ac14854..7a51ba3b8 100644
--- a/lib/atc/ATC_Coupling.h
+++ b/lib/atc/ATC_Coupling.h
@@ -1,456 +1,456 @@
 #ifndef ATC_COUPLING_H
 #define ATC_COUPLING_H
 
 #include <set>
 #include <map>
 #include <string>
 #include <utility>
 
 // ATC headers
 #include "ATC_Method.h"
 #include "ExtrinsicModel.h"
 
 namespace ATC {
 
   // Forward declarations
   class PrescribedDataManager;
   class AtomicRegulator;
   class TimeIntegrator;
   class ReferencePositions;
 
   /**
    *  @class  ATC_Coupling
    *  @brief  Base class for atom-continuum coupling 
    */
 
   class ATC_Coupling : public ATC_Method {
 
   public: /** methods */
     
     friend class ExtrinsicModel; // friend is not inherited
     friend class ExtrinsicModelTwoTemperature;
     friend class ExtrinsicModelDriftDiffusion; 
     friend class ExtrinsicModelDriftDiffusionConvection;
     friend class ExtrinsicModelElectrostatic;
     friend class ExtrinsicModelElectrostaticMomentum;
     friend class SchrodingerPoissonSolver;
     friend class SliceSchrodingerPoissonSolver;
 
     /** constructor */
     ATC_Coupling(std::string groupName, double **& perAtomArray, LAMMPS_NS::Fix * thisFix);
 
     /** destructor */
     virtual ~ATC_Coupling();
 
     /** parser/modifier */
     virtual bool modify(int narg, char **arg);
 
     /** pre neighbor */
     virtual void pre_neighbor();
 
     /** pre exchange */
     virtual void pre_exchange();
     virtual void reset_atoms(){};
 
     /** pre force */
     virtual void pre_force();
 
     /** post force */
     virtual void post_force();
 
     /** pre integration run */
     virtual void initialize();
 
     /** flags whether a methods reset is required */
     virtual bool reset_methods() const;
 
     /** post integration run : called at end of run or simulation */
     virtual void finish();
 
     /** first time, before atomic integration */
     virtual void pre_init_integrate();
 
     /** Predictor phase, Verlet first step for velocity and position */
     virtual void init_integrate(); 
 
     /** Predictor phase, executed after Verlet */
     virtual void post_init_integrate();
 
     /** Corrector phase, executed after Verlet*/
     
     virtual void post_final_integrate();
 
     /** pre/post atomic force calculation in minimize */
     virtual void min_pre_force(){};
     virtual void min_post_force(){};
 
     // data access
     /** get map general atomic shape function matrix to overlap region */
     SPAR_MAT &get_atom_to_overlap_mat() {return atomToOverlapMat_.set_quantity();};
     /** get map general atomic shape function matrix to overlap region */
     SPAR_MAN &atom_to_overlap_mat() {return atomToOverlapMat_;};
     /** check if atomic quadrature is being used for MD_ONLY nodes */
     bool atom_quadrature_on(){return atomQuadForInternal_;};
     const std::set<std::string> & boundary_face_names() {return boundaryFaceNames_;};
     /** access to boundary integration method */
     int boundary_integration_type() {return bndyIntType_;};
     void set_boundary_integration_type(int boundaryIntegrationType) 
     {bndyIntType_ = boundaryIntegrationType;};
     void set_boundary_face_set(const std::set< std::pair<int,int> > * boundaryFaceSet) 
     {bndyFaceSet_ = boundaryFaceSet;};
     BoundaryIntegrationType parse_boundary_integration
       (int narg, char **arg, const std::set< std::pair<int,int> > * boundaryFaceSet);
     TemperatureDefType temperature_def() const {return temperatureDef_;};
     void set_temperature_def(TemperatureDefType tdef) {temperatureDef_ = tdef;};
 
 //--------------------------------------------------------
 
     /** access to all boundary fluxes */
     FIELDS &boundary_fluxes() {return boundaryFlux_;}; 
     /** wrapper for FE_Engine's compute_boundary_flux functions */
     void compute_boundary_flux(const Array2D<bool> & rhs_mask,
                                const FIELDS &fields, 
                                FIELDS &rhs,
                                const Array< std::set <int> > atomMaterialGroups,
                                const VectorDependencyManager<SPAR_MAT * > * shpFcnDerivs,
                                const SPAR_MAN * shpFcn = NULL,
                                const DIAG_MAN * atomicWeights = NULL,
                                const MatrixDependencyManager<DenseMatrix, bool> * elementMask = NULL,
                                const RegulatedNodes * nodeSet = NULL);
     /** access to full right hand side / forcing vector */
     FIELDS &rhs() {return rhs_;};
 
     DENS_MAN &field_rhs(FieldName thisField) { return rhs_[thisField]; };
     /** allow FE_Engine to construct ATC structures after mesh is constructed */
     virtual void initialize_mesh_data(void);  
 
 // public for FieldIntegrator
     bool source_atomic_quadrature(FieldName field)  
-      { return (sourceIntegration_ == FULL_DOMAIN_ATOMIC_QUADRATURE_SOURCE); }
+      { return false; }
     ATC::IntegrationDomainType source_integration() 
       { return sourceIntegration_; }
 
     /** wrapper for FE_Engine's compute_sources */
     void compute_sources_at_atoms(const RHS_MASK & rhsMask,
                                   const FIELDS & fields,
                                   const PhysicsModel * physicsModel,
                                   FIELD_MATS & atomicSources);
     /** computes tangent matrix using atomic quadrature near FE region */
    void masked_atom_domain_rhs_tangent(const std::pair<FieldName,FieldName> row_col,
                                        const RHS_MASK & rhsMask,      
                                        const FIELDS & fields,                  
                                        SPAR_MAT & stiffness,
                                        const PhysicsModel * physicsModel);
     /** wrapper for FE_Engine's compute_rhs_vector functions */
     void compute_rhs_vector(const RHS_MASK & rhs_mask,
                             const FIELDS &fields, 
                             FIELDS &rhs,
                             const IntegrationDomainType domain, // = FULL_DOMAIN
                             const PhysicsModel * physicsModel=NULL);
    /** wrapper for FE_Engine's compute_tangent_matrix */
    void compute_rhs_tangent(const std::pair<FieldName,FieldName> row_col,
                             const RHS_MASK & rhsMask,      
                             const FIELDS & fields,                  
                             SPAR_MAT & stiffness,
                             const IntegrationDomainType integrationType,
                             const PhysicsModel * physicsModel=NULL);
 
    /** PDE type */
    WeakEquation::PDE_Type pde_type(const FieldName fieldName) const;
    /** is dynamic PDE */
    bool is_dynamic(const FieldName fieldName) const;
 
 // public for ImplicitSolveOperator
     /** return pointer to PrescribedDataManager */
     PrescribedDataManager * prescribed_data_manager() 
       { return prescribedDataMgr_; }
 // public for Kinetostat
     DIAG_MAT &get_mass_mat(FieldName thisField)
       { return massMats_[thisField].set_quantity();};
     /** */
     DENS_MAN &atomic_source(FieldName thisField){return atomicSources_[thisField];};
 
 
     //---------------------------------------------------------------
     /** \name materials  */
     //---------------------------------------------------------------
     /*@{*/
     /** access to element to material map */
     Array<int> &element_to_material_map(void){return elementToMaterialMap_;}
     /*@}*/
 
     /** check if method is tracking charge */
     bool track_charge() {return trackCharge_;};
 
     void set_mass_mat_time_filter(FieldName thisField,TimeFilterManager::FilterIntegrationType filterIntegrationType);
 
     /** return referece to ExtrinsicModelManager */
     ExtrinsicModelManager & extrinsic_model_manager() 
       { return extrinsicModelManager_; }
     /** access to time integrator */
     const TimeIntegrator * time_integrator(const FieldName & field) const {
       _ctiIt_ = timeIntegrators_.find(field);
       if (_ctiIt_ == timeIntegrators_.end()) return NULL;
       return _ctiIt_->second;
     };
 
     //---------------------------------------------------------------
     /** \name managers */
     //---------------------------------------------------------------
     /*@{*/
     /** allow FE_Engine to construct data manager after mesh is constructed */
     void construct_prescribed_data_manager (void); 
     /** method to create physics model */
     void create_physics_model(const PhysicsType & physicsType,
                               std::string matFileName);
     /** access to physics model */
     PhysicsModel * physics_model() {return physicsModel_; };
     /*@}*/
 
     //---------------------------------------------------------------
     /** \name creation */
     //---------------------------------------------------------------
     /*@{*/
     /** set up atom to material identification */
     virtual void reset_atom_materials();
     /** */
     void reset_node_mask();
     /** */
     void reset_overlap_map();
     /*@}*/
 
     //---------------------------------------------------------------
     /** \name output/restart */
     //---------------------------------------------------------------
     /*@{*/
     void pack_fields(RESTART_LIST & data);
     virtual void  read_restart_data(std::string fileName_, RESTART_LIST & data);
     virtual void write_restart_data(std::string fileName_, RESTART_LIST & data);
     void output() { ATC_Method::output(); }
     /*@}*/
 
     //---------------------------------------------------------------
     /** \name initial & boundary conditions  */
     //---------------------------------------------------------------
     /*@{*/
     /** mask for computation of fluxes */
     void set_fixed_nodes();
     /** set initial conditions by changing fields */
     void set_initial_conditions();
     /*@}*/
 
     //---------------------------------------------------------------
     /** \name sources  */
     //---------------------------------------------------------------
     /** calculate and set matrix of sources_ */
     void set_sources();
     /** assemble various contributions to the heat flux in the atomic region */
     void compute_atomic_sources(const Array2D<bool> & rhs_mask,
                                 const FIELDS &fields, 
                                 FIELDS &atomicSources);
 
     DENS_MAT &get_source(FieldName thisField){return sources_[thisField].set_quantity();};
     DENS_MAN &source(FieldName thisField){return sources_[thisField];};
     FIELDS & sources(){return sources_;};
     /** access to name atomic source terms */
     DENS_MAT &get_atomic_source(FieldName thisField){return atomicSources_[thisField].set_quantity();};
     /** access to name extrinsic source terms */
     DENS_MAT &get_extrinsic_source(FieldName thisField){return extrinsicSources_[thisField].set_quantity();};
     DENS_MAN &extrinsic_source(FieldName thisField){return extrinsicSources_[thisField];};
 
     /** nodal projection of a field through the physics model */
 
     void nodal_projection(const FieldName & fieldName,
                           const PhysicsModel * physicsModel,
                           FIELD & field);
     /*@}*/
 
     //---------------------------------------------------------------
     /** \name fluxes  */
     //---------------------------------------------------------------
     /*@{*/
     /** access for field mask */
     Array2D<bool> &field_mask() {return fieldMask_;}; 
     /** create field mask */
     void reset_flux_mask();
     /** field mask for intrinsic integration */
     Array2D<bool> intrinsicMask_;
     /** wrapper for FE_Engine's compute_flux functions */
     void compute_flux(const Array2D<bool> & rhs_mask,
                       const FIELDS &fields, 
                       GRAD_FIELD_MATS &flux,
                       const PhysicsModel * physicsModel=NULL);
     /** evaluate rhs on the atomic domain which is near the FE region */
     void masked_atom_domain_rhs_integral(const Array2D<bool> & rhs_mask,
                                          const FIELDS &fields, 
                                          FIELDS &rhs,
                                          const PhysicsModel * physicsModel);
     /** evaluate rhs on a specified domain defined by mask and physics model */
     void evaluate_rhs_integral(const Array2D<bool> & rhs_mask,
                            const FIELDS &fields, 
                            FIELDS &rhs,
                            const IntegrationDomainType domain,
                            const PhysicsModel * physicsModel=NULL);
     /** access to boundary fluxes */
     DENS_MAT &get_boundary_flux(FieldName thisField){return boundaryFlux_[thisField].set_quantity();};
     DENS_MAN &boundary_flux(FieldName thisField){return boundaryFlux_[thisField];};
     /** access to finite element right-hand side data */
     DENS_MAT &get_field_rhs(FieldName thisField)
       { return rhs_[thisField].set_quantity(); };
     /*@}*/
 
     //---------------------------------------------------------------
     /** \name mass matrices  */
     //---------------------------------------------------------------
     /*@{*/
     // atomic field time derivative filtering
     virtual void init_filter(void);
     // mass matrix filtering
     void delete_mass_mat_time_filter(FieldName thisField);
     /** compute mass matrix for requested field */
     void compute_mass_matrix(FieldName thisField, PhysicsModel * physicsModel = NULL);
     /** updates filtering of MD contributions */
     void update_mass_matrix(FieldName thisField);
     /** compute the mass matrix components coming from MD integration */
     virtual void compute_md_mass_matrix(FieldName thisField,
                                         DIAG_MAT & massMats);
 
   private: /** methods */
     ATC_Coupling(); // do not define
 
   protected: /** data */
 
     //---------------------------------------------------------------
     /** initialization routines */
     //---------------------------------------------------------------
     /** sets up all data necessary to define the computational geometry */
     virtual void set_computational_geometry();
     /** constructs all data which is updated with time integration, i.e. fields */
     virtual void construct_time_integration_data();
     /** create methods, e.g. time integrators, filters */
     virtual void construct_methods();
     /** set up data which is dependency managed */
     virtual void construct_transfers();
     /** sets up mol transfers */
     virtual void construct_molecule_transfers();
     /** sets up accumulant & interpolant */
     virtual void construct_interpolant();
     /** reset number of local atoms */
     virtual void reset_nlocal();
 
     //---------------------------------------------------------------
     /** status */
     //---------------------------------------------------------------
     /*@{*/
     /** flag on if FE nodes in MD region should be initialized to projected MD values */
     bool consistentInitialization_;
     bool equilibriumStart_;
     bool useFeMdMassMatrix_;
     /** flag to determine if charge is tracked */
     bool trackCharge_;
     /** temperature definition model */
     TemperatureDefType temperatureDef_;
     /*@}*/
 
     //---------------------------------------------------------------
     /** \name managers */
     //---------------------------------------------------------------
     /*@{*/
     /** prescribed data handler */
     PrescribedDataManager * prescribedDataMgr_;
     /** pointer to physics model */
     PhysicsModel * physicsModel_;
     /** manager for extrinsic models */
     ExtrinsicModelManager extrinsicModelManager_;
     /** manager for regulator */
     AtomicRegulator * atomicRegulator_;
     /** managers for time integrators per field */
     std::map<FieldName,TimeIntegrator * > timeIntegrators_;
     /** time integrator iterator */
     mutable std::map<FieldName,TimeIntegrator * >::iterator _tiIt_;
     /** time integrator const iterator */
     mutable std::map<FieldName,TimeIntegrator * >::const_iterator _ctiIt_;
     /*@}*/
 
     //---------------------------------------------------------------
     /** materials */
     //---------------------------------------------------------------
     /*@{*/
     Array<int> elementToMaterialMap_;  // ATOMIC_TAG * elementToMaterialMap_;
     /** atomic ATC material tag */
     
     
     Array< std::set <int> > atomMaterialGroups_;  // ATOMIC_TAG*atomMaterialGroups_;
     Array< std::set <int> > atomMaterialGroupsMask_;  // ATOMIC_TAG*atomMaterialGroupsMask_;
     /*@}*/
 
     //---------------------------------------------------------------
     /** computational geometry */
     //---------------------------------------------------------------
     /*@{*/
     bool atomQuadForInternal_;  
     MatrixDependencyManager<DenseMatrix, bool> * elementMask_;
     MatrixDependencyManager<DenseMatrix, bool> * elementMaskMass_;
     MatrixDependencyManager<DenseMatrix, bool> * elementMaskMassMd_;
     /** operator to compute the mass matrix for the momentum equation from MD integration */
     AtfShapeFunctionRestriction * nodalAtomicMass_;
     /** operator to compute the dimensionless mass matrix from MD integration */
     AtfShapeFunctionRestriction * nodalAtomicCount_;
     /** operator to compute mass matrix from MD */
     AtfShapeFunctionRestriction * nodalAtomicHeatCapacity_;
     MatrixDependencyManager<DenseMatrix, bool> * create_full_element_mask();
     MatrixDependencyManager<DenseMatrix, int> * create_element_set_mask(const std::string & elementSetName);
     LargeToSmallAtomMap * internalToMask_;
     MatrixDependencyManager<DenseMatrix, int> * internalElement_;
     MatrixDependencyManager<DenseMatrix, int> * ghostElement_;
     DenseMatrixTransfer<int> * nodalGeometryType_;
     /*@}*/
 
     /** \name boundary integration */
     /*@{*/
     /** boundary flux quadrature */
     int bndyIntType_;
     const std::set< std::pair<int,int> > * bndyFaceSet_;
     std::set<std::string> boundaryFaceNames_;
     /*@}*/
 
     //----------------------------------------------------------------
     /** \name shape function matrices */
     //----------------------------------------------------------------
     /*@{*/
     DIAG_MAN * atomicWeightsMask_;
     SPAR_MAN * shpFcnMask_;
     VectorDependencyManager<SPAR_MAT * > * shpFcnDerivsMask_;
     Array<bool> atomMask_;  
     SPAR_MAN atomToOverlapMat_;
     DIAG_MAN nodalMaskMat_;
     /*@}*/
 
     //---------------------------------------------------------------
     /** \name PDE data */
     //---------------------------------------------------------------
     /*@{*/
     /** mask for computation of fluxes */
     Array2D<bool> fieldMask_; 
     DIAG_MAT fluxMask_;
     DIAG_MAT fluxMaskComplement_;
     /** sources */
     FIELDS sources_;
     FIELDS atomicSources_;
     FIELDS extrinsicSources_;
     ATC::IntegrationDomainType sourceIntegration_;
     SPAR_MAT stiffnessAtomDomain_;
     /** rhs/forcing terms */
     FIELDS rhs_;  // for pde
     FIELDS rhsAtomDomain_; // for thermostat
     FIELDS boundaryFlux_; // for thermostat & rhs pde
     // DATA structures for tracking individual species and molecules
     FIELD_POINTERS atomicFields_;
     /*@}*/
 
     // workspace variables
     mutable DENS_MAT _deltaQuantity_;
   };
 };
 
 #endif
diff --git a/lib/atc/ATC_TypeDefs.h b/lib/atc/ATC_TypeDefs.h
index ab24a8a22..8906d5c63 100644
--- a/lib/atc/ATC_TypeDefs.h
+++ b/lib/atc/ATC_TypeDefs.h
@@ -1,626 +1,625 @@
 #ifndef ATC_TYPEDEFS_H
 #define ATC_TYPEDEFS_H
 
 #include <set>
 #include <vector>
 #include <map>
 #include <utility>
 #include <string>
 
 #ifdef NEW_LAMMPS
 #include "lmptype.h"
 #endif
 
 #include "Array.h"
 #include "Array2D.h"
 
 #include "MatrixLibrary.h"
 #include "DependencyManager.h"
 
 namespace ATC
 {
   /** physical constants */
   static const double kBeV_ = 8.617343e-5;// [eV/K]
 
   /** unsigned ints, when needed */
   typedef int INDEX; 
 
   /** elementset integral */
   enum ElementsetOperationType {
     ELEMENTSET_TOTAL=0,
     ELEMENTSET_AVERAGE
   };
   /** faceset integral */
   enum FacesetIntegralType {
     BOUNDARY_INTEGRAL=0,
     CONTOUR_INTEGRAL
   };
 
   /** nodeset operation */
   enum NodesetOperationType {
     NODESET_SUM=0,
     NODESET_AVERAGE
   };
 
   /** boundary integration */
   enum BoundaryIntegrationType {
     NO_QUADRATURE=0,
     FE_QUADRATURE,
     FE_INTERPOLATION
   };
   /** domain integration */
   enum IntegrationDomainType {
     FULL_DOMAIN=0,
     ATOM_DOMAIN,
-    FE_DOMAIN,
-    FULL_DOMAIN_ATOMIC_QUADRATURE_SOURCE 
+    FE_DOMAIN
   };
   /** domain decomposition */
   enum DomainDecompositionType {
     REPLICATED_MEMORY=0,
     DISTRIBUTED_MEMORY 
   };
   /** atomic weight specification */
   enum AtomicWeightType {
     USER=0,
     LATTICE,
     ELEMENT,
     REGION,
     GROUP,
     MULTISCALE,
     NODE,
     NODE_ELEMENT,
     READ_IN
   };
   /** geometry location with respect to MD domain */
   enum GeometryType {
     FE_ONLY = 0,
     MD_ONLY,
     BOUNDARY
   };
   /** enumerated type for atomic reference frame */
   enum AtomToElementMapType {
     LAGRANGIAN=0,
     EULERIAN
   };
   /* enumerated type for coupling matrix structure */
   enum MatrixStructure {
     FULL=0,     // contributions from all nodes
     LOCALIZED,  // contributions only from nodes with sources
     LUMPED      // row-sum lumped version of full matrix
   };
   /* enumerated type for distinguishing ghost from internal atoms */
   enum AtomType {
     INTERNAL=0,
     GHOST,
     ALL,
     PROC_GHOST,
     NO_ATOMS,
     NUM_ATOM_TYPES
   };
   /** field types */
   enum FieldName { 
       TIME=-2,
       POSITION=-1,
       TEMPERATURE=0, // Intrinsic Fields
       DISPLACEMENT,
       VELOCITY,
       MASS_DENSITY,
       CHARGE_DENSITY,
       SPECIES_CONCENTRATION,
       ELECTRON_DENSITY, // Extrinsic Fields
       ELECTRON_VELOCITY,
       ELECTRON_TEMPERATURE,
       ELECTRIC_POTENTIAL,
       ELECTRON_WAVEFUNCTION,
       ELECTRON_WAVEFUNCTIONS, 
       ELECTRON_WAVEFUNCTION_ENERGIES, 
       FERMI_ENERGY, 
       MOMENTUM,
       PROJECTED_VELOCITY,
       KINETIC_TEMPERATURE,
       THERMAL_ENERGY,
       KINETIC_ENERGY,
       STRESS,
       KINETIC_STRESS,
       HEAT_FLUX,
       CHARGE_FLUX,
       SPECIES_FLUX,
       INTERNAL_ENERGY,
       REFERENCE_POTENTIAL_ENERGY,
       POTENTIAL_ENERGY,
       ENERGY,
       NUMBER_DENSITY,
       ESHELBY_STRESS,
       CAUCHY_BORN_STRESS,
       CAUCHY_BORN_ENERGY,
       CAUCHY_BORN_ESHELBY_STRESS,
       TRANSFORMED_STRESS,
       VACANCY_CONCENTRATION,
       ROTATION,
       STRETCH,
       DIPOLE_MOMENT,
       QUADRUPOLE_MOMENT,
       CAUCHY_BORN_ELASTIC_DEFORMATION_GRADIENT,
       DISLOCATION_DENSITY,
       NUM_TOTAL_FIELDS 
   };
   const int NUM_FIELDS = ELECTRON_WAVEFUNCTION+1; 
 
 #define NDIM 3
   static const int FieldSizes[NUM_TOTAL_FIELDS] = {
     1, // TEMPERATURE
     NDIM, // DISPLACEMENT
     NDIM, // VELOCITY
     1, // MASS_DENSITY
     1, // CHARGE_DENSITY
     0, // SPECIES_CONCENTRATION - VARIABLE
     1, // ELECTRON_DENSITY
     NDIM, // ELECTRON_VELOCITY
     1, // ELECTRON_TEMPERATURE
     1, // ELECTRIC_POTENTIAL
     1, // ELECTRON_WAVEFUNCTION ?
     0, // ELECTRON_WAVEFUNCTIONS - VARIABLE
     0, // ELECTRON_WAVEFUNCTION_ENERGIES - VARIABLE
     1, // FERMI_ENERGY
     NDIM, // MOMENTUM
     NDIM, // PROJECTED_VELOCITY
     1, // KINETIC_TEMPERATURE
     1, // THERMAL_ENERGY
     1, // KINETIC_ENERGY
     NDIM*NDIM, // STRESS
     NDIM*NDIM, // KINETIC_STRESS
     NDIM, // HEAT_FLUX
     NDIM, // CHARGE_FLUX
     0, // SPECIES_FLUX - VARIABLE
     1, // INTERNAL_ENERGY
     1, // REFERENCE_POTENTIAL_ENERGY
     1, // POTENTIAL_ENERGY
     1, // ENERGY
     1, // NUMBER_DENSITY
     NDIM*NDIM, // ESHELBY_STRESS
     NDIM*NDIM, // CAUCHY_BORN_STRESS,
     1, //  CAUCHY_BORN_ENERGY,
     NDIM*NDIM, //  CAUCHY_BORN_ESHELBY_STRESS,
     NDIM*NDIM, //  TRANSFORMED_STRESS,
     1, //  VACANCY_CONCENTRATION,
     NDIM*NDIM, //  ROTATION,
     NDIM*NDIM, //  STRETCH,
     NDIM, // DIPOLE_MOMENT,
     NDIM, // QUADRUPOLE_MOMENT,
     NDIM*NDIM, //  CAUCHY_BORN_ELASTIC_DEFORMATION_GRADIENT,
     NDIM*NDIM //  DISLOCATION_DENSITY
   };
 
   enum NodalAtomicFieldNormalization { 
     NO_NORMALIZATION=0,
     VOLUME_NORMALIZATION, NUMBER_NORMALIZATION, MASS_NORMALIZATION,
     MASS_MATRIX
   };
 
   inline FieldName use_mass_matrix(FieldName in) {
     if (in == TEMPERATURE) return in;
     else                   return MASS_DENSITY;
   }
 
   /** enums for FE Element and Interpolate classes */
   enum FeEltGeometry   {HEXA, TETRA};
   enum FeIntQuadrature {NODAL, GAUSS1, GAUSS2, GAUSS3, FACE};
 
   /** field name enum to string */
   inline FeIntQuadrature string_to_FIQ(const std::string &str) 
   {
     if (str == "nodal")
       return NODAL;
     else if (str == "gauss1")
       return GAUSS1;
     else if (str == "gauss2")
       return GAUSS2;
     else if (str == "gauss3")
       return GAUSS3;
     else if (str == "face")
       return FACE;
     else
       throw ATC_Error("Bad quadrature input" + str + ".");
   }
 
   /** field name enum to string */
   inline std::string field_to_string(const FieldName index) 
   {
     switch (index) {
     case TEMPERATURE:
       return "temperature";
     case DISPLACEMENT:
       return "displacement";
     case VELOCITY:
       return "velocity";
     case MASS_DENSITY:
       return "mass_density";
     case CHARGE_DENSITY:
       return "charge_density";
     case ELECTRON_DENSITY:
       return "electron_density";
     case ELECTRON_VELOCITY:
       return "electron_velocity";
     case ELECTRON_TEMPERATURE:
       return "electron_temperature";
     case ELECTRIC_POTENTIAL:
       return "electric_potential";
     case ELECTRON_WAVEFUNCTION:
       return "electron_wavefunction";
     case ELECTRON_WAVEFUNCTIONS:
       return "electron_wavefunctions";
     case ELECTRON_WAVEFUNCTION_ENERGIES:
       return "electron_wavefunction_energies";
     case FERMI_ENERGY:
       return "fermi_energy";
     case MOMENTUM:
       return "momentum";
     case PROJECTED_VELOCITY:
       return "projected_velocity";
     case KINETIC_TEMPERATURE:
       return "kinetic_temperature";
     case THERMAL_ENERGY:
       return "thermal_energy";
     case KINETIC_ENERGY:
       return "kinetic_energy";
     case STRESS:
       return "stress";
     case KINETIC_STRESS:
       return "kinetic_stress";
     case ESHELBY_STRESS:
       return "eshelby_stress";
     case CAUCHY_BORN_STRESS:
       return "cauchy_born_stress";
     case CAUCHY_BORN_ENERGY:
       return "cauchy_born_energy";
     case CAUCHY_BORN_ESHELBY_STRESS:
       return "cauchy_born_eshelby_stress";
     case HEAT_FLUX:
       return "heat_flux";
     case CHARGE_FLUX:
       return "charge_flux";
     case SPECIES_FLUX:
       return "species_flux";
     case INTERNAL_ENERGY:
       return "internal_energy";
     case POTENTIAL_ENERGY:
       return "potential_energy";
     case REFERENCE_POTENTIAL_ENERGY:
       return "reference_potential_energy";
     case ENERGY:
       return "energy";
     case NUMBER_DENSITY:
       return "number_density";
     case TRANSFORMED_STRESS:
       return "transformed_stress";
     case VACANCY_CONCENTRATION:
       return "vacancy_concentration";
     case SPECIES_CONCENTRATION:
       return "species_concentration";
     case ROTATION:
       return "rotation";
     case STRETCH:
       return "stretch";
     case DIPOLE_MOMENT:
       return "dipole_moment";
     case QUADRUPOLE_MOMENT:
       return "quadrupole_moment";
     default:
       throw ATC_Error("field not found in field_to_string");
     }
   };
 
   /** string to field enum */
   inline FieldName string_to_field(const std::string & name) 
   {
     if      (name=="temperature")
       return TEMPERATURE;
     else if (name=="displacement")
       return DISPLACEMENT;
     else if (name=="velocity")
       return VELOCITY;
     else if (name=="mass_density")
       return MASS_DENSITY;
     else if (name=="charge_density")
       return CHARGE_DENSITY;
     else if (name=="electron_density")
       return ELECTRON_DENSITY;
     else if (name=="electron_velocity")
       return ELECTRON_VELOCITY;
     else if (name=="electron_temperature")
       return ELECTRON_TEMPERATURE;
     else if (name=="electric_potential")
       return ELECTRIC_POTENTIAL;
     else if (name=="electron_wavefunction")
       return ELECTRON_WAVEFUNCTION;
     else if (name=="electron_wavefunctions")
       return ELECTRON_WAVEFUNCTIONS;
     else if (name=="electron_wavefunction_energies")
       return ELECTRON_WAVEFUNCTION_ENERGIES;
     else if (name=="fermi_energy")
       return FERMI_ENERGY;
     else if (name=="momentum")
       return MOMENTUM;
     else if (name=="projected_velocity")
       return PROJECTED_VELOCITY;
     else if (name=="kinetic_temperature")
       return KINETIC_TEMPERATURE; // temperature from total KE
     else if (name=="thermal_energy")
       return THERMAL_ENERGY;
     else if (name=="kinetic_energy")
       return KINETIC_ENERGY;
     else if (name=="stress")
       return STRESS;
     else if (name=="kinetic_stress")
       return KINETIC_STRESS;
     else if (name=="eshelby_stress")
       return ESHELBY_STRESS;
     else if (name=="cauchy_born_stress")
       return CAUCHY_BORN_STRESS;
     else if (name=="cauchy_born_energy")
       return CAUCHY_BORN_ENERGY;
     else if (name=="cauchy_born_eshelby_stress")
       return CAUCHY_BORN_ESHELBY_STRESS;
     else if (name=="heat_flux")
       return HEAT_FLUX;
     else if (name=="charge_flux")
       return CHARGE_FLUX;
     else if (name=="species_flux")
       return SPECIES_FLUX;
     else if (name=="internal_energy")
       return INTERNAL_ENERGY;
     else if (name=="reference_potential_energy")
       return REFERENCE_POTENTIAL_ENERGY;
     else if (name=="potential_energy")
       return POTENTIAL_ENERGY;
     else if (name=="energy")
       return ENERGY;
     else if (name=="number_density")
       return NUMBER_DENSITY;
     else if (name=="transformed_stress")
       return TRANSFORMED_STRESS;
     else if (name=="vacancy_concentration")
       return VACANCY_CONCENTRATION;
     else if (name=="species_concentration")
       return SPECIES_CONCENTRATION;
     else if (name=="rotation")
       return ROTATION;
     else if (name=="stretch")
       return STRETCH;
     else if (name=="dipole_moment")
       return DIPOLE_MOMENT;
     else if (name=="quadrupole_moment")
       return QUADRUPOLE_MOMENT;
     else
       throw ATC_Error(name + " is not a valid field");
   };
 
   inline bool is_intrinsic(const FieldName & field_enum) 
   {
     if  (field_enum==TEMPERATURE
       || field_enum==DISPLACEMENT
       || field_enum==VELOCITY
       || field_enum==MASS_DENSITY
       || field_enum==CHARGE_DENSITY
       || field_enum==SPECIES_CONCENTRATION
       || field_enum==KINETIC_TEMPERATURE
       || field_enum==POTENTIAL_ENERGY
       || field_enum==REFERENCE_POTENTIAL_ENERGY
      )   return true;
     else return false;
   };
 
   inline std::string field_to_intrinsic_name(const FieldName index) 
   {
     if (is_intrinsic(index)) {
       return "NodalAtomic"+ATC_Utility::to_cap(field_to_string(index));
     }
     else {
       throw ATC_Error("field "+field_to_string(index)+" is not an intrinsic field");
     }
   }
   inline std::string field_to_restriction_name(const FieldName index) 
   {
     if (is_intrinsic(index)) {
       return "Restricted"+ATC_Utility::to_cap(field_to_string(index));
     }
     else {
       throw ATC_Error("field "+field_to_string(index)+" is not an intrinsic field");
     }
   }
   inline std::string field_to_prolongation_name(const FieldName index) 
   {
     return "Prolonged"+ATC_Utility::to_cap(field_to_string(index));
   }
 
 
   /** types of temperature definitions */
   enum TemperatureDefType {
     NONE = 0,
     KINETIC,
     TOTAL
   };
 
   /** string to temperature definition enum */
   inline bool string_to_temperature_def(const std::string & name, TemperatureDefType & index) {
     if      (name=="none")
       index = NONE;
     else if (name=="kinetic")
       index = KINETIC;
     else if (name=="total")
       index = TOTAL;
     else {
       throw ATC_Error("temperature operator type "+name+" not valid");
       return false;
     }
 
     return true;
   };
 
   /** solver types */
   enum SolverType { DIRECT=0, ITERATIVE};
   enum DirichletType {DIRICHLET_PENALTY=0, DIRICHLET_CONDENSE};
 
   /** physics types */
   enum PhysicsType
   {
     NO_PHYSICS=0, // for post-processing only
     THERMAL,
     ELASTIC,
     SHEAR,
     THERMO_ELASTIC,
     SPECIES // aka Mass
   };
  
   /** rhs types */
   enum FluxType 
   {
     FLUX = 0,           // has a source weighted by gradient of shape function
     SOURCE,             // has a source term weighted by the shape function
     PRESCRIBED_SOURCE,  // has a prescribed source term
     ROBIN_SOURCE,       // has a Robin source term
     EXTRINSIC_SOURCE,   // has an extrinsic source term
     NUM_FLUX
   };
 
   /** stiffness/ derivative of rhs types */
   enum StiffnessType 
   {
     BB_STIFFNESS = 0,    
     NN_STIFFNESS,    
     BN_STIFFNESS,    
     NB_STIFFNESS,    
     NUM_STIFFNESS
   };
 
   /** LAMMPS atom type identifiers */
   enum IdType {
     ATOM_TYPE=0,
     ATOM_GROUP
   };
 
   /** molecule size identifiers */
   enum MolSize {
     MOL_SMALL=0,
     MOL_LARGE
   };
 
   /** basic */
   typedef std::pair<int, int> PAIR;
 
   /** typedefs for compact set of bc values */
   typedef std::set < std::pair < int, double> > BC_SET; // node, value
   typedef std::vector< BC_SET > BCS;  // dof (node, value)
   typedef std::set<int> NSET;  // nodeset
   typedef std::set<PAIR> FSET;  // faceset
   typedef std::set<int> ESET;  // elementset
 
   /** typedefs for N and B integrand functions */
   typedef std::set<FieldName> ARG_NAMES; 
   typedef std::map<FieldName, ATC_matrix::DenseMatrix<double> > ARGS; 
   typedef ATC::MatrixDependencyManager<ATC_matrix::DenseMatrix, double>  FIELD;
   typedef std::vector<ATC::MatrixDependencyManager<ATC_matrix::DenseMatrix, double> >  GRAD_FIELD;
   typedef std::map<FieldName, ATC::MatrixDependencyManager<ATC_matrix::DenseMatrix, double> > FIELDS;
   typedef std::map<FieldName, ATC::MatrixDependencyManager<ATC_matrix::DenseMatrix, double> * > FIELD_POINTERS;
   typedef std::map<FieldName, ATC_matrix::DenseMatrix<double> > FIELD_MATS;
   typedef std::map<std::string, ATC::MatrixDependencyManager<ATC_matrix::DenseMatrix, double> > TAG_FIELDS;
   typedef std::map<FieldName, std::vector<ATC::MatrixDependencyManager<ATC_matrix::DenseMatrix, double> > > GRAD_FIELDS;
   typedef std::map<FieldName, std::vector<ATC_matrix::DenseMatrix<double> > > GRAD_FIELD_MATS;
   typedef std::map<FieldName, ATC::MatrixDependencyManager<DiagonalMatrix, double> > MASS_MATS;
   typedef std::map<FieldName, ATC::MatrixDependencyManager<SparseMatrix, double> > CON_MASS_MATS;
   typedef ATC::MatrixDependencyManager<ATC_matrix::DenseMatrix, double> DENS_MAN;
   typedef ATC::MatrixDependencyManager<SparseMatrix, double> SPAR_MAN;
   typedef ATC::MatrixDependencyManager<ParSparseMatrix, double> PAR_SPAR_MAN;
   typedef ATC::MatrixDependencyManager<DiagonalMatrix, double> DIAG_MAN;
   typedef ATC::MatrixDependencyManager<ParDiagonalMatrix, double> PAR_DIAG_MAN;
 
   /** typedefs for WeakEquation evaluation */
   typedef Array2D<bool>  RHS_MASK;
 
   /** typedefs for input/output */
   typedef std::map<std::string, const ATC_matrix::Matrix<double>*> OUTPUT_LIST;
   typedef std::map<std::string, ATC_matrix::Matrix<double>*> RESTART_LIST;
 
   typedef  std::pair<int, int>  ID_PAIR;
   typedef std::vector< std::pair<int, int> > ID_LIST;
 
   /** misc typedefs */
   class XT_Function;
   class UXT_Function;
   typedef std::map<FieldName, std::map<PAIR, Array<XT_Function*> > > SURFACE_SOURCE;
   typedef std::map<FieldName, std::map<PAIR, Array<UXT_Function*> > > ROBIN_SURFACE_SOURCE;
   typedef std::map<FieldName, Array2D<XT_Function *> > VOLUME_SOURCE;
   typedef std::map<std::string, ATC::MatrixDependencyManager<ATC_matrix::DenseMatrix, double> > ATOMIC_DATA;
   
   /** typedefs for FE_Mesh */
   typedef std::map<std::string, std::set<int > > NODE_SET_MAP;
   typedef std::map<std::string, std::set<int > > ELEMENT_SET_MAP;
   typedef std::map<std::string, std::set<PAIR> > FACE_SET_MAP;
 
   /** string to index */
   // inline vs. static is used to avoid compiler warnings that the function isn't used
   // the compiler seems to just check if static functions are used in the file they're
   // declared in rather than all the files that include the header,
   // same for arrays (but not primitives, e.g. ints) hopefully this also speeds up the code
   inline bool string_to_index(const std::string & dim, int & index, int & sgn)
   {
     char dir;
     if (dim.empty()) return false;
     sgn = (dim[0] == '-') ? -1 : 1;
     dir = dim[dim.size()-1]; // dir is last character
     if      (dir == 'x') index = 0;
     else if (dir == 'y') index = 1;
     else if (dir == 'z') index = 2;
     else return false;
     return true;
   };
 
   /** string to index */
   inline std::string index_to_string(const int &index)
   {
     if      (index==0) return "x";
     else if (index==1) return "y";
     else if (index==2) return "z";
     return "unknown";
   };
 
   /** string to index */
   inline bool string_to_index(const std::string &dim, int &index)
   {
     if (dim=="x")
       index = 0;
     else if (dim=="y")
       index = 1;
     else if (dim=="z")
       index = 2;
     else
       return false;
 
     return true;
   };
 
   inline std::string print_mask(const Array2D<bool> & rhsMask) 
   {
     std::string msg;
     for (int i = 0; i < NUM_FIELDS; i++) {       
       FieldName field = (FieldName) i;
       std::string name = field_to_string(field);
       if (rhsMask(field,FLUX)
           || rhsMask(field,SOURCE)           
           || rhsMask(field,PRESCRIBED_SOURCE)
           || rhsMask(field,ROBIN_SOURCE)
           || rhsMask(field,EXTRINSIC_SOURCE))  {
         msg = "RHS_MASK: " + name;
         if (rhsMask(field,FLUX))              msg += " flux";
         if (rhsMask(field,SOURCE))            msg += " source";
         if (rhsMask(field,PRESCRIBED_SOURCE)) msg += " prescribed_src";
         if (rhsMask(field,ROBIN_SOURCE))      msg += " robin_src";
         if (rhsMask(field,EXTRINSIC_SOURCE))  msg += " extrinsic_src";
       }
     }
     return msg;
   }
 }
 
 #endif
diff --git a/lib/atc/PoissonSolver.cpp b/lib/atc/PoissonSolver.cpp
index c2c062ede..a5b4668e5 100644
--- a/lib/atc/PoissonSolver.cpp
+++ b/lib/atc/PoissonSolver.cpp
@@ -1,216 +1,213 @@
 #include "PoissonSolver.h"
 #include "ATC_Coupling.h"
 #include "FE_Engine.h"
 #include "PhysicsModel.h"
 #include "PrescribedDataManager.h"
 #include "LinearSolver.h"
 #include <utility>
 #include <iostream>
 
 using std::pair;
 
 
 
 namespace ATC {
 
 // ====================================================================
 //  PoissonSolver
 // ====================================================================
 PoissonSolver::PoissonSolver(
     const FieldName fieldName,
     const PhysicsModel * physicsModel,
     const FE_Engine * feEngine,
     const PrescribedDataManager * prescribedDataMgr,
     /*const*/ ATC_Coupling * atc,
     const Array2D<bool> & rhsMask,
     const int solverType,
     bool parallel
 )
   : atc_(atc),
     feEngine_(feEngine),
     prescribedDataMgr_(prescribedDataMgr),
     physicsModel_(physicsModel),
     fieldName_(fieldName),
     rhsMask_(rhsMask),
     linear_(false),
     solver_(NULL),
     solverNL_(NULL),
     tangent_(NULL),
     solverType_(solverType),
     solverTol_(0),
     solverMaxIter_(0),
     integrationType_(FULL_DOMAIN),
     parallel_(parallel)
 {
   if (physicsModel_->has_linear_rhs(fieldName)) {
     linear_ = true;
     rhsMask_(fieldName,FLUX) = false; 
   }
 
   else {
     rhsMask_(fieldName,FLUX)   = true; 
     rhsMask_(fieldName,SOURCE) = true; 
   }
   if (prescribedDataMgr_->has_robin_source(fieldName)) {
     
     
     
     rhsMask_(fieldName,ROBIN_SOURCE) = true;
   }
 }
 // --------------------------------------------------------------------
 PoissonSolver::~PoissonSolver() 
 { 
   if (tangent_) delete tangent_;
   if (solverNL_) delete solverNL_;
   if (solver_) delete solver_;
 }
 
 // --------------------------------------------------------------------
 //  Parser
 // --------------------------------------------------------------------
 
 bool PoissonSolver::modify(int narg, char **arg)
 {
   bool match = false;
   /*! \page man_poisson_solver fix_modify AtC poisson_solver 
       \section syntax
       fix_modify AtC poisson_solver mesh create <nx> <ny> <nz> <region-id> 
       <f|p> <f|p> <f|p>
       - nx ny nz = number of elements in x, y, z
       - region-id = id of region that is to be meshed
       - f p p  = perioidicity flags for x, y, z
       \section examples
       <TT> fix_modify AtC poisson_solver mesh create 10 1 1 feRegion p p p </TT>
       \section description
       Creates a uniform mesh in a rectangular region
       \section restrictions
       creates only uniform rectangular grids in a rectangular region
       \section related
       \section default
       none
   */
   int argIdx = 0;
   if (strcmp(arg[argIdx],"poisson_solver")==0) {
     argIdx++;
     if (strcmp(arg[argIdx],"mesh")==0) {
       argIdx++;
       // create a FE_Engine
 
       //feEngine_ = new FE_Engine(this); need alternate constructor?
       // send args to new engine
       // arg[0] = "mesh";
       // arg[1] = "create";
       // feEngine_->modify(narg,arg);
 
     }
   } 
   return match;
 }
 // --------------------------------------------------------------------
 //  Initialize
 // --------------------------------------------------------------------
 void PoissonSolver::initialize(void)
 {
   nNodes_ = feEngine_->num_nodes();
 
-  if (atc_->source_atomic_quadrature(fieldName_))  
-    integrationType_ = FULL_DOMAIN_ATOMIC_QUADRATURE_SOURCE;
-
   // compute penalty for Dirichlet boundaries
   if (prescribedDataMgr_->none_fixed(fieldName_))  
     throw ATC_Error("Poisson solver needs Dirichlet data");
 
   const BC_SET & bcs = (prescribedDataMgr_->bcs(fieldName_))[0];
 
   if (linear_) { // constant rhs 
     if (! solver_ ) {
       pair<FieldName,FieldName> row_col(fieldName_,fieldName_);
       Array2D <bool> rhsMask(NUM_FIELDS,NUM_FLUX);
       rhsMask = false; rhsMask(fieldName_,FLUX) = true;
       if (prescribedDataMgr_->has_robin_source(fieldName_)) {
         rhsMask(fieldName_,ROBIN_SOURCE) = true;
       }
       // compute stiffness for Poisson solve
       atc_->compute_rhs_tangent(row_col, rhsMask, atc_->fields(),
         stiffness_, FULL_DOMAIN, physicsModel_);
       // create solver
       solver_ = new LinearSolver(stiffness_,bcs,solverType_,-1,parallel_);
     }
     else {
 
     // re-initialize
     solver_->initialize(&bcs);
     }
     if (solverTol_) solver_->set_tolerance(solverTol_);
     if (solverMaxIter_) solver_->set_max_iterations(solverMaxIter_);
 
   }
   else {
 //  print_mask(rhsMask_);
     if ( solverNL_ )  delete solverNL_;
     tangent_ = new PhysicsModelTangentOperator(atc_,physicsModel_, rhsMask_, integrationType_, fieldName_); 
 
     solverNL_ = new NonLinearSolver(tangent_,&bcs,0,parallel_);
     
     if (solverTol_) solverNL_->set_residual_tolerance(solverTol_);
     if (solverMaxIter_) solverNL_->set_max_iterations(solverMaxIter_);
   }
 }
 
 // --------------------------------------------------------------------
 //  Solve
 // --------------------------------------------------------------------
 bool PoissonSolver::solve(FIELDS & fields, FIELDS & rhs)
 {
   atc_->compute_rhs_vector(rhsMask_, fields, rhs,
     integrationType_, physicsModel_);
   CLON_VEC f = column(fields[fieldName_].set_quantity(),0);
   CLON_VEC r = column(rhs[fieldName_].quantity(),0);
   bool converged = false;
   if (linear_) {converged = solver_->solve(f,r);}
   else         {converged = solverNL_->solve(f);}
 
   if (atc_->source_atomic_quadrature(fieldName_) 
     && LammpsInterface::instance()->atom_charge() ) set_charges(fields);
   return converged;
 }
 bool PoissonSolver::solve(DENS_MAT & field, const DENS_MAT & rhs) 
 {
 
   CLON_VEC f = column(field,0);
   CLON_VEC r = column(rhs,0);
   bool converged = false;
   if (linear_) {converged = solver_->solve(f,r);}
   else         {converged = solverNL_->solve(f);}
 
   if (atc_->source_atomic_quadrature(fieldName_) 
     && LammpsInterface::instance()->atom_charge() ) set_charges(atc_->fields());
   return converged;
 }
 
 // --------------------------------------------------------------------
 //  set charges on atoms
 // --------------------------------------------------------------------
 void PoissonSolver::set_charges(FIELDS & fields)
 {
   FIELD_MATS sources;
   
   atc_->compute_sources_at_atoms(rhsMask_, fields, physicsModel_,sources);
   FIELD_MATS::const_iterator nField = sources.find(fieldName_);
   if (nField != sources.end()) {
     const DENS_MAT & electronCharges = nField->second;
     double *  q = LammpsInterface::instance()->atom_charge();
     int nLocal = atc_->nlocal();
     if (nLocal > 0) {
       const Array<int> & i2a = atc_->internal_to_atom_map();
       for (int i=0; i < nLocal; i++) {
 
         int atomIdx = i2a(i);
         q[atomIdx] = -electronCharges(i,0);
       }
     }
   }
 }
 
 } // namespace ATC