diff --git a/src/Makefile.shlib b/src/Makefile.shlib
index 17c1bdb8f..a4770269d 100644
--- a/src/Makefile.shlib
+++ b/src/Makefile.shlib
@@ -1,43 +1,43 @@
 # LAMMPS shared library multiple-machine Makefile
 
 SHELL = /bin/sh
 
 # Definitions
 
 ROOT =	lammps
 EXE =	lib$(ROOT)_$@.so
 
-SRC =	angle.cpp angle_charmm.cpp angle_class2.cpp angle_cosine.cpp angle_cosine_delta.cpp angle_cosine_periodic.cpp angle_cosine_squared.cpp angle_harmonic.cpp angle_hybrid.cpp angle_table.cpp atom.cpp atom_map.cpp atom_vec.cpp atom_vec_angle.cpp atom_vec_atomic.cpp atom_vec_bond.cpp atom_vec_charge.cpp atom_vec_dipole.cpp atom_vec_ellipsoid.cpp atom_vec_full.cpp atom_vec_hybrid.cpp atom_vec_line.cpp atom_vec_molecular.cpp atom_vec_peri.cpp atom_vec_sphere.cpp atom_vec_tri.cpp balance.cpp bond.cpp bond_class2.cpp bond_fene.cpp bond_fene_expand.cpp bond_harmonic.cpp bond_hybrid.cpp bond_morse.cpp bond_nonlinear.cpp bond_quartic.cpp bond_table.cpp change_box.cpp comm.cpp compute.cpp compute_angle_local.cpp compute_atom_molecule.cpp compute_bond_local.cpp compute_centro_atom.cpp compute_cluster_atom.cpp compute_cna_atom.cpp compute_com.cpp compute_com_molecule.cpp compute_contact_atom.cpp compute_coord_atom.cpp compute_damage_atom.cpp compute_dihedral_local.cpp compute_displace_atom.cpp compute_erotate_asphere.cpp compute_erotate_sphere.cpp compute_erotate_sphere_atom.cpp compute_event_displace.cpp compute_group_group.cpp compute_gyration.cpp compute_gyration_molecule.cpp compute_heat_flux.cpp compute_improper_local.cpp compute_inertia_molecule.cpp compute_ke.cpp compute_ke_atom.cpp compute_msd.cpp compute_msd_molecule.cpp compute_pair.cpp compute_pair_local.cpp compute_pe.cpp compute_pe_atom.cpp compute_pressure.cpp compute_property_atom.cpp compute_property_local.cpp compute_property_molecule.cpp compute_rdf.cpp compute_reduce.cpp compute_reduce_region.cpp compute_slice.cpp compute_stress_atom.cpp compute_temp.cpp compute_temp_asphere.cpp compute_temp_com.cpp compute_temp_deform.cpp compute_temp_partial.cpp compute_temp_profile.cpp compute_temp_ramp.cpp compute_temp_region.cpp compute_temp_sphere.cpp compute_ti.cpp create_atoms.cpp create_box.cpp delete_atoms.cpp delete_bonds.cpp dihedral.cpp dihedral_charmm.cpp dihedral_class2.cpp dihedral_harmonic.cpp dihedral_helix.cpp dihedral_hybrid.cpp dihedral_multi_harmonic.cpp dihedral_opls.cpp displace_atoms.cpp domain.cpp dump.cpp dump_atom.cpp dump_cfg.cpp dump_custom.cpp dump_dcd.cpp dump_image.cpp dump_local.cpp dump_xtc.cpp dump_xyz.cpp error.cpp finish.cpp fix.cpp fix_adapt.cpp fix_addforce.cpp fix_append_atoms.cpp fix_ave_atom.cpp fix_ave_correlate.cpp fix_ave_histo.cpp fix_ave_spatial.cpp fix_ave_time.cpp fix_aveforce.cpp fix_balance.cpp fix_bond_break.cpp fix_bond_create.cpp fix_bond_swap.cpp fix_box_relax.cpp fix_deform.cpp fix_deposit.cpp fix_drag.cpp fix_dt_reset.cpp fix_efield.cpp fix_enforce2d.cpp fix_evaporate.cpp fix_event.cpp fix_event_prd.cpp fix_event_tad.cpp fix_external.cpp fix_freeze.cpp fix_gcmc.cpp fix_gravity.cpp fix_heat.cpp fix_indent.cpp fix_langevin.cpp fix_lineforce.cpp fix_minimize.cpp fix_momentum.cpp fix_move.cpp fix_msst.cpp fix_neb.cpp fix_nh.cpp fix_nh_asphere.cpp fix_nh_sphere.cpp fix_nph.cpp fix_nph_asphere.cpp fix_nph_sphere.cpp fix_nphug.cpp fix_npt.cpp fix_npt_asphere.cpp fix_npt_sphere.cpp fix_nve.cpp fix_nve_asphere.cpp fix_nve_asphere_noforce.cpp fix_nve_limit.cpp fix_nve_line.cpp fix_nve_noforce.cpp fix_nve_sphere.cpp fix_nve_tri.cpp fix_nvt.cpp fix_nvt_asphere.cpp fix_nvt_sllod.cpp fix_nvt_sphere.cpp fix_orient_fcc.cpp fix_peri_neigh.cpp fix_planeforce.cpp fix_pour.cpp fix_press_berendsen.cpp fix_print.cpp fix_qeq_comb.cpp fix_read_restart.cpp fix_recenter.cpp fix_respa.cpp fix_restrain.cpp fix_rigid.cpp fix_rigid_nh.cpp fix_rigid_nph.cpp fix_rigid_npt.cpp fix_rigid_nve.cpp fix_rigid_nvt.cpp fix_setforce.cpp fix_shake.cpp fix_shear_history.cpp fix_spring.cpp fix_spring_rg.cpp fix_spring_self.cpp fix_srd.cpp fix_store_force.cpp fix_store_state.cpp fix_temp_berendsen.cpp fix_temp_rescale.cpp fix_thermal_conductivity.cpp fix_tmd.cpp fix_ttm.cpp fix_viscosity.cpp fix_viscous.cpp fix_wall.cpp fix_wall_colloid.cpp fix_wall_gran.cpp fix_wall_harmonic.cpp fix_wall_lj126.cpp fix_wall_lj93.cpp fix_wall_piston.cpp fix_wall_reflect.cpp fix_wall_region.cpp fix_wall_srd.cpp force.cpp group.cpp image.cpp improper.cpp improper_class2.cpp improper_cvff.cpp improper_harmonic.cpp improper_hybrid.cpp improper_umbrella.cpp input.cpp integrate.cpp irregular.cpp kspace.cpp lammps.cpp lattice.cpp library.cpp  math_extra.cpp memory.cpp min.cpp min_cg.cpp min_fire.cpp min_hftn.cpp min_linesearch.cpp min_quickmin.cpp min_sd.cpp minimize.cpp modify.cpp neb.cpp neigh_bond.cpp neigh_derive.cpp neigh_full.cpp neigh_gran.cpp neigh_half_bin.cpp neigh_half_multi.cpp neigh_half_nsq.cpp neigh_list.cpp neigh_request.cpp neigh_respa.cpp neigh_stencil.cpp neighbor.cpp output.cpp pair.cpp pair_adp.cpp pair_airebo.cpp pair_beck.cpp pair_bop.cpp pair_born.cpp pair_born_coul_wolf.cpp pair_brownian.cpp pair_brownian_poly.cpp pair_buck.cpp pair_buck_coul_cut.cpp pair_colloid.cpp pair_comb.cpp pair_coul_cut.cpp pair_coul_debye.cpp pair_coul_dsf.cpp pair_coul_wolf.cpp pair_dipole_cut.cpp pair_dpd.cpp pair_dpd_tstat.cpp pair_dsmc.cpp pair_eam.cpp pair_eam_alloy.cpp pair_eam_alloy_opt.cpp pair_eam_fs.cpp pair_eam_fs_opt.cpp pair_eam_opt.cpp pair_eim.cpp pair_gauss.cpp pair_gayberne.cpp pair_gran_hertz_history.cpp pair_gran_hooke.cpp pair_gran_hooke_history.cpp pair_hbond_dreiding_lj.cpp pair_hbond_dreiding_morse.cpp pair_hybrid.cpp pair_hybrid_overlay.cpp pair_lcbop.cpp pair_line_lj.cpp pair_lj96_cut.cpp pair_lj_charmm_coul_charmm.cpp pair_lj_charmm_coul_charmm_implicit.cpp pair_lj_class2.cpp pair_lj_class2_coul_cut.cpp pair_lj_class2_coul_long.cpp pair_lj_cubic.cpp pair_lj_cut.cpp pair_lj_cut_coul_cut.cpp pair_lj_cut_coul_debye.cpp pair_lj_cut_coul_dsf.cpp pair_lj_cut_opt.cpp pair_lj_expand.cpp pair_lj_gromacs.cpp pair_lj_gromacs_coul_gromacs.cpp pair_lj_smooth.cpp pair_lj_smooth_linear.cpp pair_lubricate.cpp pair_lubricateU.cpp pair_lubricateU_poly.cpp pair_lubricate_poly.cpp pair_mie_cut.cpp pair_morse.cpp pair_morse_opt.cpp pair_peri_lps.cpp pair_peri_pmb.cpp pair_rebo.cpp pair_resquared.cpp pair_soft.cpp pair_sw.cpp pair_table.cpp pair_tersoff.cpp pair_tersoff_zbl.cpp pair_tri_lj.cpp pair_yukawa.cpp pair_yukawa_colloid.cpp prd.cpp procmap.cpp random_mars.cpp random_park.cpp read_data.cpp read_dump.cpp read_restart.cpp reader.cpp reader_native.cpp reader_xyz.cpp region.cpp region_block.cpp region_cone.cpp region_cylinder.cpp region_intersect.cpp region_plane.cpp region_prism.cpp region_sphere.cpp region_union.cpp replicate.cpp rerun.cpp respa.cpp run.cpp set.cpp special.cpp tad.cpp temper.cpp thermo.cpp timer.cpp universe.cpp update.cpp variable.cpp velocity.cpp verlet.cpp verlet_split.cpp write_restart.cpp xdr_compat.cpp 
+SRC =	angle.cpp angle_charmm.cpp angle_class2.cpp angle_cosine.cpp angle_cosine_delta.cpp angle_cosine_periodic.cpp angle_cosine_squared.cpp angle_harmonic.cpp angle_hybrid.cpp angle_table.cpp atom.cpp atom_map.cpp atom_vec.cpp atom_vec_angle.cpp atom_vec_atomic.cpp atom_vec_body.cpp atom_vec_bond.cpp atom_vec_charge.cpp atom_vec_dipole.cpp atom_vec_ellipsoid.cpp atom_vec_full.cpp atom_vec_hybrid.cpp atom_vec_line.cpp atom_vec_molecular.cpp atom_vec_peri.cpp atom_vec_sphere.cpp atom_vec_tri.cpp balance.cpp body.cpp body_nparticle.cpp bond.cpp bond_class2.cpp bond_fene.cpp bond_fene_expand.cpp bond_harmonic.cpp bond_hybrid.cpp bond_morse.cpp bond_nonlinear.cpp bond_quartic.cpp bond_table.cpp change_box.cpp comm.cpp compute.cpp compute_angle_local.cpp compute_atom_molecule.cpp compute_body_local.cpp compute_bond_local.cpp compute_centro_atom.cpp compute_cluster_atom.cpp compute_cna_atom.cpp compute_com.cpp compute_com_molecule.cpp compute_contact_atom.cpp compute_coord_atom.cpp compute_damage_atom.cpp compute_dihedral_local.cpp compute_displace_atom.cpp compute_erotate_asphere.cpp compute_erotate_sphere.cpp compute_erotate_sphere_atom.cpp compute_event_displace.cpp compute_group_group.cpp compute_gyration.cpp compute_gyration_molecule.cpp compute_heat_flux.cpp compute_improper_local.cpp compute_inertia_molecule.cpp compute_ke.cpp compute_ke_atom.cpp compute_msd.cpp compute_msd_molecule.cpp compute_pair.cpp compute_pair_local.cpp compute_pe.cpp compute_pe_atom.cpp compute_pressure.cpp compute_property_atom.cpp compute_property_local.cpp compute_property_molecule.cpp compute_rdf.cpp compute_reduce.cpp compute_reduce_region.cpp compute_slice.cpp compute_stress_atom.cpp compute_temp.cpp compute_temp_asphere.cpp compute_temp_com.cpp compute_temp_deform.cpp compute_temp_partial.cpp compute_temp_profile.cpp compute_temp_ramp.cpp compute_temp_region.cpp compute_temp_sphere.cpp compute_ti.cpp create_atoms.cpp create_box.cpp delete_atoms.cpp delete_bonds.cpp dihedral.cpp dihedral_charmm.cpp dihedral_class2.cpp dihedral_harmonic.cpp dihedral_helix.cpp dihedral_hybrid.cpp dihedral_multi_harmonic.cpp dihedral_opls.cpp displace_atoms.cpp domain.cpp dump.cpp dump_atom.cpp dump_cfg.cpp dump_custom.cpp dump_dcd.cpp dump_image.cpp dump_local.cpp dump_xtc.cpp dump_xyz.cpp error.cpp finish.cpp fix.cpp fix_adapt.cpp fix_addforce.cpp fix_append_atoms.cpp fix_ave_atom.cpp fix_ave_correlate.cpp fix_ave_histo.cpp fix_ave_spatial.cpp fix_ave_time.cpp fix_aveforce.cpp fix_balance.cpp fix_bond_break.cpp fix_bond_create.cpp fix_bond_swap.cpp fix_box_relax.cpp fix_deform.cpp fix_deposit.cpp fix_drag.cpp fix_dt_reset.cpp fix_efield.cpp fix_enforce2d.cpp fix_evaporate.cpp fix_event.cpp fix_event_prd.cpp fix_event_tad.cpp fix_external.cpp fix_freeze.cpp fix_gcmc.cpp fix_gravity.cpp fix_heat.cpp fix_indent.cpp fix_langevin.cpp fix_lineforce.cpp fix_minimize.cpp fix_momentum.cpp fix_move.cpp fix_msst.cpp fix_neb.cpp fix_nh.cpp fix_nh_asphere.cpp fix_nh_sphere.cpp fix_nph.cpp fix_nph_asphere.cpp fix_nph_sphere.cpp fix_nphug.cpp fix_npt.cpp fix_npt_asphere.cpp fix_npt_sphere.cpp fix_nve.cpp fix_nve_asphere.cpp fix_nve_asphere_noforce.cpp fix_nve_body.cpp fix_nve_limit.cpp fix_nve_line.cpp fix_nve_noforce.cpp fix_nve_sphere.cpp fix_nve_tri.cpp fix_nvt.cpp fix_nvt_asphere.cpp fix_nvt_sllod.cpp fix_nvt_sphere.cpp fix_orient_fcc.cpp fix_peri_neigh.cpp fix_planeforce.cpp fix_pour.cpp fix_press_berendsen.cpp fix_print.cpp fix_qeq_comb.cpp fix_read_restart.cpp fix_recenter.cpp fix_respa.cpp fix_restrain.cpp fix_rigid.cpp fix_rigid_nh.cpp fix_rigid_nph.cpp fix_rigid_npt.cpp fix_rigid_nve.cpp fix_rigid_nvt.cpp fix_rigid_small.cpp fix_setforce.cpp fix_shake.cpp fix_shear_history.cpp fix_spring.cpp fix_spring_rg.cpp fix_spring_self.cpp fix_srd.cpp fix_store_force.cpp fix_store_state.cpp fix_temp_berendsen.cpp fix_temp_rescale.cpp fix_thermal_conductivity.cpp fix_tmd.cpp fix_ttm.cpp fix_viscosity.cpp fix_viscous.cpp fix_wall.cpp fix_wall_colloid.cpp fix_wall_gran.cpp fix_wall_harmonic.cpp fix_wall_lj126.cpp fix_wall_lj93.cpp fix_wall_piston.cpp fix_wall_reflect.cpp fix_wall_region.cpp fix_wall_srd.cpp force.cpp group.cpp image.cpp improper.cpp improper_class2.cpp improper_cvff.cpp improper_harmonic.cpp improper_hybrid.cpp improper_umbrella.cpp input.cpp integrate.cpp irregular.cpp kspace.cpp lammps.cpp lattice.cpp library.cpp  math_extra.cpp memory.cpp min.cpp min_cg.cpp min_fire.cpp min_hftn.cpp min_linesearch.cpp min_quickmin.cpp min_sd.cpp minimize.cpp modify.cpp neb.cpp neigh_bond.cpp neigh_derive.cpp neigh_full.cpp neigh_gran.cpp neigh_half_bin.cpp neigh_half_multi.cpp neigh_half_nsq.cpp neigh_list.cpp neigh_request.cpp neigh_respa.cpp neigh_stencil.cpp neighbor.cpp output.cpp pair.cpp pair_adp.cpp pair_airebo.cpp pair_beck.cpp pair_body.cpp pair_bop.cpp pair_born.cpp pair_born_coul_wolf.cpp pair_brownian.cpp pair_brownian_poly.cpp pair_buck.cpp pair_buck_coul_cut.cpp pair_colloid.cpp pair_comb.cpp pair_coul_cut.cpp pair_coul_debye.cpp pair_coul_dsf.cpp pair_coul_wolf.cpp pair_dipole_cut.cpp pair_dpd.cpp pair_dpd_tstat.cpp pair_dsmc.cpp pair_eam.cpp pair_eam_alloy.cpp pair_eam_alloy_opt.cpp pair_eam_fs.cpp pair_eam_fs_opt.cpp pair_eam_opt.cpp pair_eim.cpp pair_gauss.cpp pair_gayberne.cpp pair_gran_hertz_history.cpp pair_gran_hooke.cpp pair_gran_hooke_history.cpp pair_hbond_dreiding_lj.cpp pair_hbond_dreiding_morse.cpp pair_hybrid.cpp pair_hybrid_overlay.cpp pair_lcbop.cpp pair_line_lj.cpp pair_lj96_cut.cpp pair_lj_charmm_coul_charmm.cpp pair_lj_charmm_coul_charmm_implicit.cpp pair_lj_class2.cpp pair_lj_class2_coul_cut.cpp pair_lj_class2_coul_long.cpp pair_lj_cubic.cpp pair_lj_cut.cpp pair_lj_cut_coul_cut.cpp pair_lj_cut_coul_debye.cpp pair_lj_cut_coul_dsf.cpp pair_lj_cut_opt.cpp pair_lj_cut_tip4p_cut.cpp pair_lj_expand.cpp pair_lj_gromacs.cpp pair_lj_gromacs_coul_gromacs.cpp pair_lj_smooth.cpp pair_lj_smooth_linear.cpp pair_lubricate.cpp pair_lubricateU.cpp pair_lubricateU_poly.cpp pair_lubricate_poly.cpp pair_mie_cut.cpp pair_morse.cpp pair_morse_opt.cpp pair_peri_lps.cpp pair_peri_pmb.cpp pair_rebo.cpp pair_resquared.cpp pair_soft.cpp pair_sw.cpp pair_table.cpp pair_tersoff.cpp pair_tersoff_zbl.cpp pair_tri_lj.cpp pair_yukawa.cpp pair_yukawa_colloid.cpp prd.cpp procmap.cpp random_mars.cpp random_park.cpp read_data.cpp read_dump.cpp read_restart.cpp reader.cpp reader_native.cpp reader_xyz.cpp region.cpp region_block.cpp region_cone.cpp region_cylinder.cpp region_intersect.cpp region_plane.cpp region_prism.cpp region_sphere.cpp region_union.cpp replicate.cpp rerun.cpp respa.cpp run.cpp set.cpp special.cpp tad.cpp temper.cpp thermo.cpp timer.cpp universe.cpp update.cpp variable.cpp velocity.cpp verlet.cpp verlet_split.cpp write_restart.cpp xdr_compat.cpp 
 
-INC =	accelerator_cuda.h accelerator_omp.h angle.h angle_charmm.h angle_class2.h angle_cosine.h angle_cosine_delta.h angle_cosine_periodic.h angle_cosine_squared.h angle_harmonic.h angle_hybrid.h angle_table.h atom.h atom_map.h atom_masks.h atom_vec.h atom_vec_angle.h atom_vec_atomic.h atom_vec_bond.h atom_vec_charge.h atom_vec_dipole.h atom_vec_ellipsoid.h atom_vec_full.h atom_vec_hybrid.h atom_vec_line.h atom_vec_molecular.h atom_vec_peri.h atom_vec_sphere.h atom_vec_tri.h balance.h bond.h bond_class2.h bond_fene.h bond_fene_expand.h bond_harmonic.h bond_hybrid.h bond_morse.h bond_nonlinear.h bond_quartic.h bond_table.h change_box.h comm.h compute.h compute_angle_local.h compute_atom_molecule.h compute_bond_local.h compute_centro_atom.h compute_cluster_atom.h compute_cna_atom.h compute_com.h compute_com_molecule.h compute_contact_atom.h compute_coord_atom.h compute_damage_atom.h compute_dihedral_local.h compute_displace_atom.h compute_erotate_asphere.h compute_erotate_sphere.h compute_erotate_sphere_atom.h compute_event_displace.h compute_group_group.h compute_gyration.h compute_gyration_molecule.h compute_heat_flux.h compute_improper_local.h compute_inertia_molecule.h compute_ke.h compute_ke_atom.h compute_msd.h compute_msd_molecule.h compute_pair.h compute_pair_local.h compute_pe.h compute_pe_atom.h compute_pressure.h compute_property_atom.h compute_property_local.h compute_property_molecule.h compute_rdf.h compute_reduce.h compute_reduce_region.h compute_slice.h compute_stress_atom.h compute_temp.h compute_temp_asphere.h compute_temp_com.h compute_temp_deform.h compute_temp_partial.h compute_temp_profile.h compute_temp_ramp.h compute_temp_region.h compute_temp_sphere.h compute_ti.h create_atoms.h create_box.h delete_atoms.h delete_bonds.h dihedral.h dihedral_charmm.h dihedral_class2.h dihedral_harmonic.h dihedral_helix.h dihedral_hybrid.h dihedral_multi_harmonic.h dihedral_opls.h displace_atoms.h domain.h dump.h dump_atom.h dump_cfg.h dump_custom.h dump_dcd.h dump_image.h dump_local.h dump_xtc.h dump_xyz.h error.h finish.h fix.h fix_adapt.h fix_addforce.h fix_append_atoms.h fix_ave_atom.h fix_ave_correlate.h fix_ave_histo.h fix_ave_spatial.h fix_ave_time.h fix_aveforce.h fix_balance.h fix_bond_break.h fix_bond_create.h fix_bond_swap.h fix_box_relax.h fix_deform.h fix_deposit.h fix_drag.h fix_dt_reset.h fix_efield.h fix_enforce2d.h fix_evaporate.h fix_event.h fix_event_prd.h fix_event_tad.h fix_external.h fix_freeze.h fix_gcmc.h fix_gravity.h fix_heat.h fix_indent.h fix_langevin.h fix_lineforce.h fix_minimize.h fix_momentum.h fix_move.h fix_msst.h fix_neb.h fix_nh.h fix_nh_asphere.h fix_nh_sphere.h fix_nph.h fix_nph_asphere.h fix_nph_sphere.h fix_nphug.h fix_npt.h fix_npt_asphere.h fix_npt_sphere.h fix_nve.h fix_nve_asphere.h fix_nve_asphere_noforce.h fix_nve_limit.h fix_nve_line.h fix_nve_noforce.h fix_nve_sphere.h fix_nve_tri.h fix_nvt.h fix_nvt_asphere.h fix_nvt_sllod.h fix_nvt_sphere.h fix_orient_fcc.h fix_peri_neigh.h fix_planeforce.h fix_pour.h fix_press_berendsen.h fix_print.h fix_qeq_comb.h fix_read_restart.h fix_recenter.h fix_respa.h fix_restrain.h fix_rigid.h fix_rigid_nh.h fix_rigid_nph.h fix_rigid_npt.h fix_rigid_nve.h fix_rigid_nvt.h fix_setforce.h fix_shake.h fix_shear_history.h fix_spring.h fix_spring_rg.h fix_spring_self.h fix_srd.h fix_store_force.h fix_store_state.h fix_temp_berendsen.h fix_temp_rescale.h fix_thermal_conductivity.h fix_tmd.h fix_ttm.h fix_viscosity.h fix_viscous.h fix_wall.h fix_wall_colloid.h fix_wall_gran.h fix_wall_harmonic.h fix_wall_lj126.h fix_wall_lj93.h fix_wall_piston.h fix_wall_reflect.h fix_wall_region.h fix_wall_srd.h force.h group.h image.h improper.h improper_class2.h improper_cvff.h improper_harmonic.h improper_hybrid.h improper_umbrella.h input.h integrate.h irregular.h kspace.h lammps.h lattice.h library.h lmptype.h lmpwindows.h math_const.h math_extra.h memory.h min.h min_cg.h min_fire.h min_hftn.h min_linesearch.h min_quickmin.h min_sd.h minimize.h modify.h neb.h neigh_bond.h neigh_derive.h neigh_full.h neigh_gran.h neigh_half_bin.h neigh_half_multi.h neigh_half_nsq.h neigh_list.h neigh_request.h neigh_respa.h neighbor.h output.h pack.h pair.h pair_adp.h pair_airebo.h pair_beck.h pair_bop.h pair_born.h pair_born_coul_wolf.h pair_brownian.h pair_brownian_poly.h pair_buck.h pair_buck_coul_cut.h pair_colloid.h pair_comb.h pair_coul_cut.h pair_coul_debye.h pair_coul_dsf.h pair_coul_wolf.h pair_dipole_cut.h pair_dpd.h pair_dpd_tstat.h pair_dsmc.h pair_eam.h pair_eam_alloy.h pair_eam_alloy_opt.h pair_eam_fs.h pair_eam_fs_opt.h pair_eam_opt.h pair_eim.h pair_gauss.h pair_gayberne.h pair_gran_hertz_history.h pair_gran_hooke.h pair_gran_hooke_history.h pair_hbond_dreiding_lj.h pair_hbond_dreiding_morse.h pair_hybrid.h pair_hybrid_overlay.h pair_lcbop.h pair_line_lj.h pair_lj96_cut.h pair_lj_charmm_coul_charmm.h pair_lj_charmm_coul_charmm_implicit.h pair_lj_class2.h pair_lj_class2_coul_cut.h pair_lj_class2_coul_long.h pair_lj_cubic.h pair_lj_cut.h pair_lj_cut_coul_cut.h pair_lj_cut_coul_debye.h pair_lj_cut_coul_dsf.h pair_lj_cut_opt.h pair_lj_expand.h pair_lj_gromacs.h pair_lj_gromacs_coul_gromacs.h pair_lj_smooth.h pair_lj_smooth_linear.h pair_lubricate.h pair_lubricateU.h pair_lubricateU_poly.h pair_lubricate_poly.h pair_mie_cut.h pair_morse.h pair_morse_opt.h pair_peri_lps.h pair_peri_pmb.h pair_rebo.h pair_resquared.h pair_soft.h pair_sw.h pair_table.h pair_tersoff.h pair_tersoff_zbl.h pair_tri_lj.h pair_yukawa.h pair_yukawa_colloid.h pointers.h prd.h procmap.h random_mars.h random_park.h read_data.h read_dump.h read_restart.h reader.h reader_native.h reader_xyz.h region.h region_block.h region_cone.h region_cylinder.h region_intersect.h region_plane.h region_prism.h region_sphere.h region_union.h replicate.h rerun.h respa.h run.h set.h special.h style_angle.h style_atom.h style_bond.h style_command.h style_compute.h style_dihedral.h style_dump.h style_fix.h style_improper.h style_integrate.h style_kspace.h style_minimize.h style_pair.h style_reader.h style_region.h suffix.h tad.h temper.h thermo.h timer.h universe.h update.h variable.h velocity.h verlet.h verlet_split.h version.h write_restart.h xdr_compat.h 
+INC =	accelerator_cuda.h accelerator_omp.h angle.h angle_charmm.h angle_class2.h angle_cosine.h angle_cosine_delta.h angle_cosine_periodic.h angle_cosine_squared.h angle_harmonic.h angle_hybrid.h angle_table.h atom.h atom_map.h atom_masks.h atom_vec.h atom_vec_angle.h atom_vec_atomic.h atom_vec_body.h atom_vec_bond.h atom_vec_charge.h atom_vec_dipole.h atom_vec_ellipsoid.h atom_vec_full.h atom_vec_hybrid.h atom_vec_line.h atom_vec_molecular.h atom_vec_peri.h atom_vec_sphere.h atom_vec_tri.h balance.h body.h body_nparticle.h bond.h bond_class2.h bond_fene.h bond_fene_expand.h bond_harmonic.h bond_hybrid.h bond_morse.h bond_nonlinear.h bond_quartic.h bond_table.h change_box.h comm.h compute.h compute_angle_local.h compute_atom_molecule.h compute_body_local.h compute_bond_local.h compute_centro_atom.h compute_cluster_atom.h compute_cna_atom.h compute_com.h compute_com_molecule.h compute_contact_atom.h compute_coord_atom.h compute_damage_atom.h compute_dihedral_local.h compute_displace_atom.h compute_erotate_asphere.h compute_erotate_sphere.h compute_erotate_sphere_atom.h compute_event_displace.h compute_group_group.h compute_gyration.h compute_gyration_molecule.h compute_heat_flux.h compute_improper_local.h compute_inertia_molecule.h compute_ke.h compute_ke_atom.h compute_msd.h compute_msd_molecule.h compute_pair.h compute_pair_local.h compute_pe.h compute_pe_atom.h compute_pressure.h compute_property_atom.h compute_property_local.h compute_property_molecule.h compute_rdf.h compute_reduce.h compute_reduce_region.h compute_slice.h compute_stress_atom.h compute_temp.h compute_temp_asphere.h compute_temp_com.h compute_temp_deform.h compute_temp_partial.h compute_temp_profile.h compute_temp_ramp.h compute_temp_region.h compute_temp_sphere.h compute_ti.h create_atoms.h create_box.h delete_atoms.h delete_bonds.h dihedral.h dihedral_charmm.h dihedral_class2.h dihedral_harmonic.h dihedral_helix.h dihedral_hybrid.h dihedral_multi_harmonic.h dihedral_opls.h displace_atoms.h domain.h dump.h dump_atom.h dump_cfg.h dump_custom.h dump_dcd.h dump_image.h dump_local.h dump_xtc.h dump_xyz.h error.h finish.h fix.h fix_adapt.h fix_addforce.h fix_append_atoms.h fix_ave_atom.h fix_ave_correlate.h fix_ave_histo.h fix_ave_spatial.h fix_ave_time.h fix_aveforce.h fix_balance.h fix_bond_break.h fix_bond_create.h fix_bond_swap.h fix_box_relax.h fix_deform.h fix_deposit.h fix_drag.h fix_dt_reset.h fix_efield.h fix_enforce2d.h fix_evaporate.h fix_event.h fix_event_prd.h fix_event_tad.h fix_external.h fix_freeze.h fix_gcmc.h fix_gravity.h fix_heat.h fix_indent.h fix_langevin.h fix_lineforce.h fix_minimize.h fix_momentum.h fix_move.h fix_msst.h fix_neb.h fix_nh.h fix_nh_asphere.h fix_nh_sphere.h fix_nph.h fix_nph_asphere.h fix_nph_sphere.h fix_nphug.h fix_npt.h fix_npt_asphere.h fix_npt_sphere.h fix_nve.h fix_nve_asphere.h fix_nve_asphere_noforce.h fix_nve_body.h fix_nve_limit.h fix_nve_line.h fix_nve_noforce.h fix_nve_sphere.h fix_nve_tri.h fix_nvt.h fix_nvt_asphere.h fix_nvt_sllod.h fix_nvt_sphere.h fix_orient_fcc.h fix_peri_neigh.h fix_planeforce.h fix_pour.h fix_press_berendsen.h fix_print.h fix_qeq_comb.h fix_read_restart.h fix_recenter.h fix_respa.h fix_restrain.h fix_rigid.h fix_rigid_nh.h fix_rigid_nph.h fix_rigid_npt.h fix_rigid_nve.h fix_rigid_nvt.h fix_rigid_small.h fix_setforce.h fix_shake.h fix_shear_history.h fix_spring.h fix_spring_rg.h fix_spring_self.h fix_srd.h fix_store_force.h fix_store_state.h fix_temp_berendsen.h fix_temp_rescale.h fix_thermal_conductivity.h fix_tmd.h fix_ttm.h fix_viscosity.h fix_viscous.h fix_wall.h fix_wall_colloid.h fix_wall_gran.h fix_wall_harmonic.h fix_wall_lj126.h fix_wall_lj93.h fix_wall_piston.h fix_wall_reflect.h fix_wall_region.h fix_wall_srd.h force.h group.h image.h improper.h improper_class2.h improper_cvff.h improper_harmonic.h improper_hybrid.h improper_umbrella.h input.h integrate.h irregular.h kspace.h lammps.h lattice.h library.h lmptype.h lmpwindows.h math_const.h math_extra.h math_special.h memory.h min.h min_cg.h min_fire.h min_hftn.h min_linesearch.h min_quickmin.h min_sd.h minimize.h modify.h my_pool.h neb.h neigh_bond.h neigh_derive.h neigh_full.h neigh_gran.h neigh_half_bin.h neigh_half_multi.h neigh_half_nsq.h neigh_list.h neigh_request.h neigh_respa.h neighbor.h output.h pack.h pair.h pair_adp.h pair_airebo.h pair_beck.h pair_body.h pair_bop.h pair_born.h pair_born_coul_wolf.h pair_brownian.h pair_brownian_poly.h pair_buck.h pair_buck_coul_cut.h pair_colloid.h pair_comb.h pair_coul_cut.h pair_coul_debye.h pair_coul_dsf.h pair_coul_wolf.h pair_dipole_cut.h pair_dpd.h pair_dpd_tstat.h pair_dsmc.h pair_eam.h pair_eam_alloy.h pair_eam_alloy_opt.h pair_eam_fs.h pair_eam_fs_opt.h pair_eam_opt.h pair_eim.h pair_gauss.h pair_gayberne.h pair_gran_hertz_history.h pair_gran_hooke.h pair_gran_hooke_history.h pair_hbond_dreiding_lj.h pair_hbond_dreiding_morse.h pair_hybrid.h pair_hybrid_overlay.h pair_lcbop.h pair_line_lj.h pair_lj96_cut.h pair_lj_charmm_coul_charmm.h pair_lj_charmm_coul_charmm_implicit.h pair_lj_class2.h pair_lj_class2_coul_cut.h pair_lj_class2_coul_long.h pair_lj_cubic.h pair_lj_cut.h pair_lj_cut_coul_cut.h pair_lj_cut_coul_debye.h pair_lj_cut_coul_dsf.h pair_lj_cut_opt.h pair_lj_cut_tip4p_cut.h pair_lj_expand.h pair_lj_gromacs.h pair_lj_gromacs_coul_gromacs.h pair_lj_smooth.h pair_lj_smooth_linear.h pair_lubricate.h pair_lubricateU.h pair_lubricateU_poly.h pair_lubricate_poly.h pair_mie_cut.h pair_morse.h pair_morse_opt.h pair_peri_lps.h pair_peri_pmb.h pair_rebo.h pair_resquared.h pair_soft.h pair_sw.h pair_table.h pair_tersoff.h pair_tersoff_zbl.h pair_tri_lj.h pair_yukawa.h pair_yukawa_colloid.h pointers.h prd.h procmap.h random_mars.h random_park.h read_data.h read_dump.h read_restart.h reader.h reader_native.h reader_xyz.h region.h region_block.h region_cone.h region_cylinder.h region_intersect.h region_plane.h region_prism.h region_sphere.h region_union.h replicate.h rerun.h respa.h run.h set.h special.h style_angle.h style_atom.h style_body.h style_bond.h style_command.h style_compute.h style_dihedral.h style_dump.h style_fix.h style_improper.h style_integrate.h style_kspace.h style_minimize.h style_pair.h style_reader.h style_region.h suffix.h tad.h temper.h thermo.h timer.h universe.h update.h variable.h velocity.h verlet.h verlet_split.h version.h write_restart.h xdr_compat.h 
 
 OBJ =	$(SRC:.cpp=.o)
 
 # Targets
 
 help:
 	@echo 'Type "make target" where target is one of:'
 	@echo ''
 	@files="`ls MAKE/Makefile.*`"; \
 	for file in $$files; do head -1 $$file; done
 
 clean:
 	rm -rf Obj_shlib_*
 
 .DEFAULT:
 	@test -f MAKE/Makefile.$@
 	@if [ ! -d Obj_shlib_$@ ]; then mkdir Obj_shlib_$@; fi
 	@cp -p $(SRC) $(INC) Obj_shlib_$@
 	@cp MAKE/Makefile.$@ Obj_shlib_$@/Makefile
 	@if [ ! -e Makefile.package ]; \
 	  then cp Makefile.package.empty Makefile.package; fi
 	@if [ ! -e Makefile.package.settings ]; \
 	  then cp Makefile.package.settings.empty Makefile.package.settings; fi
 	@cp Makefile.package Makefile.package.settings Obj_shlib_$@
 	@cd Obj_shlib_$@; \
 	$(MAKE) $(MFLAGS) "OBJ = $(OBJ)" \
           "INC = $(INC)" "EXE = ../$(EXE)" shlib
 	@rm -f liblammps.so
 	@ln -s $(EXE) liblammps.so
 	@if [ -d Obj_shlib_$@ ]; then cd Obj_shlib_$@; \
           rm -f $(SRC) $(INC) Makefile*; fi
diff --git a/src/input.cpp b/src/input.cpp
index 51772ffec..a4cc43595 100644
--- a/src/input.cpp
+++ b/src/input.cpp
@@ -1,1506 +1,1508 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 #include "mpi.h"
 #include "stdio.h"
 #include "stdlib.h"
 #include "string.h"
 #include "ctype.h"
 #include "unistd.h"
 #include "sys/stat.h"
 #include "input.h"
 #include "style_command.h"
 #include "universe.h"
 #include "atom.h"
 #include "atom_vec.h"
 #include "comm.h"
 #include "group.h"
 #include "domain.h"
 #include "output.h"
 #include "thermo.h"
 #include "force.h"
 #include "pair.h"
 #include "min.h"
 #include "modify.h"
 #include "compute.h"
 #include "bond.h"
 #include "angle.h"
 #include "dihedral.h"
 #include "improper.h"
 #include "kspace.h"
 #include "update.h"
 #include "neighbor.h"
 #include "special.h"
 #include "variable.h"
 #include "accelerator_cuda.h"
 #include "error.h"
 #include "memory.h"
 
 #ifdef _OPENMP
 #include "omp.h"
 #endif
 
 using namespace LAMMPS_NS;
 
 #define DELTALINE 256
 #define DELTA 4
 
 /* ---------------------------------------------------------------------- */
 
 Input::Input(LAMMPS *lmp, int argc, char **argv) : Pointers(lmp)
 {
   MPI_Comm_rank(world,&me);
 
   maxline = maxcopy = maxwork = 0;
   line = copy = work = NULL;
   narg = maxarg = 0;
   arg = NULL;
 
   echo_screen = 0;
   echo_log = 1;
 
   label_active = 0;
   labelstr = NULL;
   jump_skip = 0;
 
   if (me == 0) {
     nfile = maxfile = 1;
     infiles = (FILE **) memory->smalloc(sizeof(FILE *),"input:infiles");
     infiles[0] = infile;
   } else infiles = NULL;
 
   variable = new Variable(lmp);
 
   // process command-line args
   // check for args "-var" and "-echo"
   // caller has already checked that sufficient arguments exist
 
   int iarg = 0;
   while (iarg < argc) {
     if (strcmp(argv[iarg],"-var") == 0 || strcmp(argv[iarg],"-v") == 0) {
       int jarg = iarg+3;
       while (jarg < argc && argv[jarg][0] != '-') jarg++;
       variable->set(argv[iarg+1],jarg-iarg-2,&argv[iarg+2]);
       iarg = jarg;
     } else if (strcmp(argv[iarg],"-echo") == 0 ||
                strcmp(argv[iarg],"-e") == 0) {
       narg = 1;
       char **tmp = arg;        // trick echo() into using argv instead of arg
       arg = &argv[iarg+1];
       echo();
       arg = tmp;
       iarg += 2;
     } else iarg++;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 Input::~Input()
 {
   // don't free command and arg strings
   // they just point to other allocated memory
 
   memory->sfree(line);
   memory->sfree(copy);
   memory->sfree(work);
   if (labelstr) delete [] labelstr;
   memory->sfree(arg);
   memory->sfree(infiles);
   delete variable;
 }
 
 /* ----------------------------------------------------------------------
    process all input from infile
    infile = stdin or file if command-line arg "-in" was used
 ------------------------------------------------------------------------- */
 
 void Input::file()
 {
   int m,n;
 
   while (1) {
 
     // read a line from input script
     // n = length of line including str terminator, 0 if end of file
     // if line ends in continuation char '&', concatenate next line
 
     if (me == 0) {
       m = 0;
       while (1) {
         if (maxline-m < 2) reallocate(line,maxline,0);
         if (fgets(&line[m],maxline-m,infile) == NULL) {
           if (m) n = strlen(line) + 1;
           else n = 0;
           break;
         }
         m = strlen(line);
         if (line[m-1] != '\n') continue;
 
         m--;
         while (m >= 0 && isspace(line[m])) m--;
         if (m < 0 || line[m] != '&') {
           line[m+1] = '\0';
           n = m+2;
           break;
         }
       }
     }
 
     // bcast the line
     // if n = 0, end-of-file
     // error if label_active is set, since label wasn't encountered
     // if original input file, code is done
     // else go back to previous input file
 
     MPI_Bcast(&n,1,MPI_INT,0,world);
     if (n == 0) {
       if (label_active) error->all(FLERR,"Label wasn't found in input script");
       if (me == 0) {
         if (infile != stdin) fclose(infile);
         nfile--;
       }
       MPI_Bcast(&nfile,1,MPI_INT,0,world);
       if (nfile == 0) break;
       if (me == 0) infile = infiles[nfile-1];
       continue;
     }
 
     if (n > maxline) reallocate(line,maxline,n);
     MPI_Bcast(line,n,MPI_CHAR,0,world);
 
     // echo the command unless scanning for label
 
     if (me == 0 && label_active == 0) {
       if (echo_screen && screen) fprintf(screen,"%s\n",line);
       if (echo_log && logfile) fprintf(logfile,"%s\n",line);
     }
 
     // parse the line
     // if no command, skip to next line in input script
 
     parse();
     if (command == NULL) continue;
 
     // if scanning for label, skip command unless it's a label command
 
     if (label_active && strcmp(command,"label") != 0) continue;
 
     // execute the command
 
     if (execute_command()) {
       char str[maxline+32];
       sprintf(str,"Unknown command: %s",line);
       error->all(FLERR,str);
     }
   }
 }
 
 /* ----------------------------------------------------------------------
    process all input from filename
+   called from library interface
 ------------------------------------------------------------------------- */
 
 void Input::file(const char *filename)
 {
-  // error if another nested file still open
-  // if single open file is not stdin, close it
-  // open new filename and set infile, infiles[0]
+  // error if another nested file still open, should not be possible
+  // open new filename and set infile, infiles[0], nfile
+  // call to file() will close filename and decrement nfile
 
   if (me == 0) {
     if (nfile > 1)
-      error->one(FLERR,"Another input script is already being processed");
-    if (infile != stdin) fclose(infile);
+      error->one(FLERR,"Invalid use of library file() function");
+
     infile = fopen(filename,"r");
     if (infile == NULL) {
       char str[128];
       sprintf(str,"Cannot open input script %s",filename);
       error->one(FLERR,str);
     }
     infiles[0] = infile;
-  } else infile = NULL;
+    nfile = 1;
+  }
 
   file();
 }
 
 /* ----------------------------------------------------------------------
    copy command in single to line, parse and execute it
    return command name to caller
 ------------------------------------------------------------------------- */
 
 char *Input::one(const char *single)
 {
   int n = strlen(single) + 1;
   if (n > maxline) reallocate(line,maxline,n);
   strcpy(line,single);
 
   // echo the command unless scanning for label
 
   if (me == 0 && label_active == 0) {
     if (echo_screen && screen) fprintf(screen,"%s\n",line);
     if (echo_log && logfile) fprintf(logfile,"%s\n",line);
   }
 
   // parse the line
   // if no command, just return NULL
 
   parse();
   if (command == NULL) return NULL;
 
   // if scanning for label, skip command unless it's a label command
 
   if (label_active && strcmp(command,"label") != 0) return NULL;
 
   // execute the command and return its name
 
   if (execute_command()) {
     char str[maxline+32];
     sprintf(str,"Unknown command: %s",line);
     error->all(FLERR,str);
   }
 
   return command;
 }
 
 /* ----------------------------------------------------------------------
    parse copy of command line by inserting string terminators
    strip comment = all chars from # on
    replace all $ via variable substitution
    command = first word
    narg = # of args
    arg[] = individual args
    treat text between single/double quotes as one arg
 ------------------------------------------------------------------------- */
 
 void Input::parse()
 {
   // duplicate line into copy string to break into words
 
   int n = strlen(line) + 1;
   if (n > maxcopy) reallocate(copy,maxcopy,n);
   strcpy(copy,line);
 
   // strip any # comment by replacing it with 0
   // do not strip # inside single/double quotes
 
   char quote = '\0';
   char *ptr = copy;
   while (*ptr) {
     if (*ptr == '#' && !quote) {
       *ptr = '\0';
       break;
     }
     if (*ptr == quote) quote = '\0';
     else if (*ptr == '"' || *ptr == '\'') quote = *ptr;
     ptr++;
   }
 
   // perform $ variable substitution (print changes)
   // except if searching for a label since earlier variable may not be defined
 
   if (!label_active) substitute(copy,work,maxcopy,maxwork,1);
 
   // command = 1st arg in copy string
 
   char *next;
   command = nextword(copy,&next);
   if (command == NULL) return;
 
   // point arg[] at each subsequent arg in copy string
   // nextword() inserts string terminators into copy string to delimit args
   // nextword() treats text between single/double quotes as one arg
 
   narg = 0;
   ptr = next;
   while (ptr) {
     if (narg == maxarg) {
       maxarg += DELTA;
       arg = (char **) memory->srealloc(arg,maxarg*sizeof(char *),"input:arg");
     }
     arg[narg] = nextword(ptr,&next);
     if (!arg[narg]) break;
     narg++;
     ptr = next;
   }
 }
 
 /* ----------------------------------------------------------------------
    find next word in str
    insert 0 at end of word
    ignore leading whitespace
    treat text between single/double quotes as one arg
      matching quote must be followed by whitespace char if not end of string
      strip quotes from returned word
    return ptr to start of word
    return next = ptr after word or NULL if word ended with 0
    return NULL if no word in string
 ------------------------------------------------------------------------- */
 
 char *Input::nextword(char *str, char **next)
 {
   char *start,*stop;
 
   start = &str[strspn(str," \t\n\v\f\r")];
   if (*start == '\0') return NULL;
 
   if (*start == '"' || *start == '\'') {
     stop = strchr(&start[1],*start);
     if (!stop) error->all(FLERR,"Unbalanced quotes in input line");
     if (stop[1] && !isspace(stop[1]))
       error->all(FLERR,"Input line quote not followed by whitespace");
     start++;
   } else stop = &start[strcspn(start," \t\n\v\f\r")];
 
   if (*stop == '\0') *next = NULL;
   else *next = stop+1;
   *stop = '\0';
   return start;
 }
 
 /* ----------------------------------------------------------------------
    substitute for $ variables in str using work str2 and return it
    reallocate str/str2 to hold expanded version if necessary & reset max/max2
    print updated string if flag is set and not searching for label
    label_active will be 0 if called from external class
 ------------------------------------------------------------------------- */
 
 void Input::substitute(char *&str, char *&str2, int &max, int &max2, int flag)
 {
   // use str2 as scratch space to expand str, then copy back to str
   // reallocate str and str2 as necessary
   // do not replace $ inside single/double quotes
   // var = pts at variable name, ended by NULL
   //   if $ is followed by '{', trailing '}' becomes NULL
   //   else $x becomes x followed by NULL
   // beyond = points to text following variable
 
   int i,n,paren_count;
   char immediate[256];
   char *var,*value,*beyond;
   char quote = '\0';
   char *ptr = str;
 
   n = strlen(str) + 1;
   if (n > max2) reallocate(str2,max2,n);
   *str2 = '\0';
   char *ptr2 = str2;
 
   while (*ptr) {
     // variable substitution
 
     if (*ptr == '$' && !quote) {
 
       // value = ptr to expanded variable
       // variable name between curly braces, e.g. ${a}
 
       if (*(ptr+1) == '{') {
         var = ptr+2;
         i = 0;
 
         while (var[i] != '\0' && var[i] != '}') i++;
 
         if (var[i] == '\0') error->one(FLERR,"Invalid variable name");
         var[i] = '\0';
         beyond = ptr + strlen(var) + 3;
         value = variable->retrieve(var);
 
       // immediate variable between parenthesis, e.g. $(1/2)
 
       } else if (*(ptr+1) == '(') {
         var = ptr+2;
         paren_count = 0;
         i = 0;
 
         while (var[i] != '\0' && !(var[i] == ')' && paren_count == 0)) {
           switch (var[i]) {
             case '(': paren_count++; break;
             case ')': paren_count--; break;
             default: ;
           }
           i++;
         }
 
         if (var[i] == '\0') error->one(FLERR,"Invalid immediate variable");
         var[i] = '\0';
         beyond = ptr + strlen(var) + 3;
         sprintf(immediate,"%.20g",variable->compute_equal(var));
         value = immediate;
 
       // single character variable name, e.g. $a
 
       } else {
         var = ptr;
         var[0] = var[1];
         var[1] = '\0';
         beyond = ptr + 2;
         value = variable->retrieve(var);
       }
 
       if (value == NULL) error->one(FLERR,"Substitution for illegal variable");
 
       // check if storage in str2 needs to be expanded
       // re-initialize ptr and ptr2 to the point beyond the variable.
 
       n = strlen(str2) + strlen(value) + strlen(beyond) + 1;
       if (n > max2) reallocate(str2,max2,n);
       strcat(str2,value);
       ptr2 = str2 + strlen(str2);
       ptr = beyond;
 
       // output substitution progress if requested
 
       if (flag && me == 0 && label_active == 0) {
         if (echo_screen && screen) fprintf(screen,"%s%s\n",str2,beyond);
         if (echo_log && logfile) fprintf(logfile,"%s%s\n",str2,beyond);
       }
 
       continue;
     }
 
     if (*ptr == quote) quote = '\0';
     else if (*ptr == '"' || *ptr == '\'') quote = *ptr;
 
     // copy current character into str2
 
     *ptr2++ = *ptr++;
     *ptr2 = '\0';
   }
 
   // set length of input str to length of work str2
   // copy work string back to input str
 
   if (max2 > max) reallocate(str,max,max2);
   strcpy(str,str2);
 }
 
 /* ----------------------------------------------------------------------
    rellocate a string
    if n > 0: set max >= n in increments of DELTALINE
    if n = 0: just increment max by DELTALINE
 ------------------------------------------------------------------------- */
 
 void Input::reallocate(char *&str, int &max, int n)
 {
   if (n) {
     while (n > max) max += DELTALINE;
   } else max += DELTALINE;
 
   str = (char *) memory->srealloc(str,max*sizeof(char),"input:str");
 }
 
 /* ----------------------------------------------------------------------
    process a single parsed command
    return 0 if successful, -1 if did not recognize command
 ------------------------------------------------------------------------- */
 
 int Input::execute_command()
 {
   int flag = 1;
 
   if (!strcmp(command,"clear")) clear();
   else if (!strcmp(command,"echo")) echo();
   else if (!strcmp(command,"if")) ifthenelse();
   else if (!strcmp(command,"include")) include();
   else if (!strcmp(command,"jump")) jump();
   else if (!strcmp(command,"label")) label();
   else if (!strcmp(command,"log")) log();
   else if (!strcmp(command,"next")) next_command();
   else if (!strcmp(command,"partition")) partition();
   else if (!strcmp(command,"print")) print();
   else if (!strcmp(command,"quit")) quit();
   else if (!strcmp(command,"shell")) shell();
   else if (!strcmp(command,"variable")) variable_command();
 
   else if (!strcmp(command,"angle_coeff")) angle_coeff();
   else if (!strcmp(command,"angle_style")) angle_style();
   else if (!strcmp(command,"atom_modify")) atom_modify();
   else if (!strcmp(command,"atom_style")) atom_style();
   else if (!strcmp(command,"bond_coeff")) bond_coeff();
   else if (!strcmp(command,"bond_style")) bond_style();
   else if (!strcmp(command,"boundary")) boundary();
   else if (!strcmp(command,"box")) box();
   else if (!strcmp(command,"communicate")) communicate();
   else if (!strcmp(command,"compute")) compute();
   else if (!strcmp(command,"compute_modify")) compute_modify();
   else if (!strcmp(command,"dielectric")) dielectric();
   else if (!strcmp(command,"dihedral_coeff")) dihedral_coeff();
   else if (!strcmp(command,"dihedral_style")) dihedral_style();
   else if (!strcmp(command,"dimension")) dimension();
   else if (!strcmp(command,"dump")) dump();
   else if (!strcmp(command,"dump_modify")) dump_modify();
   else if (!strcmp(command,"fix")) fix();
   else if (!strcmp(command,"fix_modify")) fix_modify();
   else if (!strcmp(command,"group")) group_command();
   else if (!strcmp(command,"improper_coeff")) improper_coeff();
   else if (!strcmp(command,"improper_style")) improper_style();
   else if (!strcmp(command,"kspace_modify")) kspace_modify();
   else if (!strcmp(command,"kspace_style")) kspace_style();
   else if (!strcmp(command,"lattice")) lattice();
   else if (!strcmp(command,"mass")) mass();
   else if (!strcmp(command,"min_modify")) min_modify();
   else if (!strcmp(command,"min_style")) min_style();
   else if (!strcmp(command,"neigh_modify")) neigh_modify();
   else if (!strcmp(command,"neighbor")) neighbor_command();
   else if (!strcmp(command,"newton")) newton();
   else if (!strcmp(command,"package")) package();
   else if (!strcmp(command,"pair_coeff")) pair_coeff();
   else if (!strcmp(command,"pair_modify")) pair_modify();
   else if (!strcmp(command,"pair_style")) pair_style();
   else if (!strcmp(command,"pair_write")) pair_write();
   else if (!strcmp(command,"processors")) processors();
   else if (!strcmp(command,"region")) region();
   else if (!strcmp(command,"reset_timestep")) reset_timestep();
   else if (!strcmp(command,"restart")) restart();
   else if (!strcmp(command,"run_style")) run_style();
   else if (!strcmp(command,"special_bonds")) special_bonds();
   else if (!strcmp(command,"suffix")) suffix();
   else if (!strcmp(command,"thermo")) thermo();
   else if (!strcmp(command,"thermo_modify")) thermo_modify();
   else if (!strcmp(command,"thermo_style")) thermo_style();
   else if (!strcmp(command,"timestep")) timestep();
   else if (!strcmp(command,"uncompute")) uncompute();
   else if (!strcmp(command,"undump")) undump();
   else if (!strcmp(command,"unfix")) unfix();
   else if (!strcmp(command,"units")) units();
 
   else flag = 0;
 
   // return if command was listed above
 
   if (flag) return 0;
 
   // check if command is added via style.h
 
   if (0) return 0;      // dummy line to enable else-if macro expansion
 
 #define COMMAND_CLASS
 #define CommandStyle(key,Class)         \
   else if (strcmp(command,#key) == 0) { \
     Class key(lmp);                     \
     key.command(narg,arg);              \
     return 0;                           \
   }
 #include "style_command.h"
 #undef COMMAND_CLASS
 
   // unrecognized command
 
   return -1;
 }
 
 /* ---------------------------------------------------------------------- */
 /* ---------------------------------------------------------------------- */
 /* ---------------------------------------------------------------------- */
 
 /* ---------------------------------------------------------------------- */
 
 void Input::clear()
 {
   if (narg > 0) error->all(FLERR,"Illegal clear command");
   lmp->destroy();
   lmp->create();
   lmp->post_create();
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::echo()
 {
   if (narg != 1) error->all(FLERR,"Illegal echo command");
 
   if (strcmp(arg[0],"none") == 0) {
     echo_screen = 0;
     echo_log = 0;
   } else if (strcmp(arg[0],"screen") == 0) {
     echo_screen = 1;
     echo_log = 0;
   } else if (strcmp(arg[0],"log") == 0) {
     echo_screen = 0;
     echo_log = 1;
   } else if (strcmp(arg[0],"both") == 0) {
     echo_screen = 1;
     echo_log = 1;
   } else error->all(FLERR,"Illegal echo command");
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::ifthenelse()
 {
   if (narg < 3) error->all(FLERR,"Illegal if command");
 
   // substitute for variables in Boolean expression for "if"
   // in case expression was enclosed in quotes
   // must substitute on copy of arg else will step on subsequent args
 
   int n = strlen(arg[0]) + 1;
   if (n > maxline) reallocate(line,maxline,n);
   strcpy(line,arg[0]);
   substitute(line,work,maxline,maxwork,0);
 
   // evaluate Boolean expression for "if"
 
   double btest = variable->evaluate_boolean(line);
 
   // bound "then" commands
 
   if (strcmp(arg[1],"then") != 0) error->all(FLERR,"Illegal if command");
 
   int first = 2;
   int iarg = first;
   while (iarg < narg &&
          (strcmp(arg[iarg],"elif") != 0 && strcmp(arg[iarg],"else") != 0))
     iarg++;
   int last = iarg-1;
 
   // execute "then" commands
   // make copies of all arg string commands
   // required because re-parsing a command via one() will wipe out args
 
   if (btest != 0.0) {
     int ncommands = last-first + 1;
     if (ncommands <= 0) error->all(FLERR,"Illegal if command");
 
     char **commands = new char*[ncommands];
     ncommands = 0;
     for (int i = first; i <= last; i++) {
       int n = strlen(arg[i]) + 1;
       if (n == 1) error->all(FLERR,"Illegal if command");
       commands[ncommands] = new char[n];
       strcpy(commands[ncommands],arg[i]);
       ncommands++;
     }
 
     for (int i = 0; i < ncommands; i++) one(commands[i]);
 
     for (int i = 0; i < ncommands; i++) delete [] commands[i];
     delete [] commands;
 
     return;
   }
 
   // done if no "elif" or "else"
 
   if (iarg == narg) return;
 
   // check "elif" or "else" until find commands to execute
   // substitute for variables and evaluate Boolean expression for "elif"
   // must substitute on copy of arg else will step on subsequent args
   // bound and execute "elif" or "else" commands
 
   while (1) {
     if (iarg+2 > narg) error->all(FLERR,"Illegal if command");
     if (strcmp(arg[iarg],"elif") == 0) {
       n = strlen(arg[iarg+1]) + 1;
       if (n > maxline) reallocate(line,maxline,n);
       strcpy(line,arg[iarg+1]);
       substitute(line,work,maxline,maxwork,0);
       btest = variable->evaluate_boolean(line);
       first = iarg+2;
     } else {
       btest = 1.0;
       first = iarg+1;
     }
 
     iarg = first;
     while (iarg < narg &&
            (strcmp(arg[iarg],"elif") != 0 && strcmp(arg[iarg],"else") != 0))
       iarg++;
     last = iarg-1;
 
     if (btest == 0.0) continue;
 
     int ncommands = last-first + 1;
     if (ncommands <= 0) error->all(FLERR,"Illegal if command");
 
     char **commands = new char*[ncommands];
     ncommands = 0;
     for (int i = first; i <= last; i++) {
       int n = strlen(arg[i]) + 1;
       if (n == 1) error->all(FLERR,"Illegal if command");
       commands[ncommands] = new char[n];
       strcpy(commands[ncommands],arg[i]);
       ncommands++;
     }
 
     // execute the list of commands
 
     for (int i = 0; i < ncommands; i++) one(commands[i]);
 
     // clean up
 
     for (int i = 0; i < ncommands; i++) delete [] commands[i];
     delete [] commands;
 
     return;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::include()
 {
   if (narg != 1) error->all(FLERR,"Illegal include command");
 
   if (me == 0) {
     if (nfile == maxfile) {
       maxfile++;
       infiles = (FILE **)
         memory->srealloc(infiles,maxfile*sizeof(FILE *),"input:infiles");
     }
     infile = fopen(arg[0],"r");
     if (infile == NULL) {
       char str[128];
       sprintf(str,"Cannot open input script %s",arg[0]);
       error->one(FLERR,str);
     }
     infiles[nfile++] = infile;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::jump()
 {
   if (narg < 1 || narg > 2) error->all(FLERR,"Illegal jump command");
 
   if (jump_skip) {
     jump_skip = 0;
     return;
   }
 
   if (me == 0) {
     if (strcmp(arg[0],"SELF") == 0) rewind(infile);
     else {
       if (infile != stdin) fclose(infile);
       infile = fopen(arg[0],"r");
       if (infile == NULL) {
         char str[128];
         sprintf(str,"Cannot open input script %s",arg[0]);
         error->one(FLERR,str);
       }
       infiles[nfile-1] = infile;
     }
   }
 
   if (narg == 2) {
     label_active = 1;
     if (labelstr) delete [] labelstr;
     int n = strlen(arg[1]) + 1;
     labelstr = new char[n];
     strcpy(labelstr,arg[1]);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::label()
 {
   if (narg != 1) error->all(FLERR,"Illegal label command");
   if (label_active && strcmp(labelstr,arg[0]) == 0) label_active = 0;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::log()
 {
   if (narg != 1) error->all(FLERR,"Illegal log command");
 
   if (me == 0) {
     if (logfile) fclose(logfile);
     if (strcmp(arg[0],"none") == 0) logfile = NULL;
     else {
       logfile = fopen(arg[0],"w");
       if (logfile == NULL) {
         char str[128];
         sprintf(str,"Cannot open logfile %s",arg[0]);
         error->one(FLERR,str);
       }
     }
     if (universe->nworlds == 1) universe->ulogfile = logfile;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::next_command()
 {
   if (variable->next(narg,arg)) jump_skip = 1;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::partition()
 {
   if (narg < 3) error->all(FLERR,"Illegal partition command");
 
   int yesflag;
   if (strcmp(arg[0],"yes") == 0) yesflag = 1;
   else if (strcmp(arg[0],"no") == 0) yesflag = 0;
   else error->all(FLERR,"Illegal partition command");
 
   int ilo,ihi;
   force->bounds(arg[1],universe->nworlds,ilo,ihi);
 
   // copy original line to copy, since will use strtok() on it
   // ptr = start of 4th word
 
   strcpy(copy,line);
   char *ptr = strtok(copy," \t\n\r\f");
   ptr = strtok(NULL," \t\n\r\f");
   ptr = strtok(NULL," \t\n\r\f");
   ptr += strlen(ptr) + 1;
   ptr += strspn(ptr," \t\n\r\f");
 
   // execute the remaining command line on requested partitions
 
   if (yesflag) {
     if (universe->iworld+1 >= ilo && universe->iworld+1 <= ihi) one(ptr);
   } else {
     if (universe->iworld+1 < ilo || universe->iworld+1 > ihi) one(ptr);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::print()
 {
   if (narg != 1) error->all(FLERR,"Illegal print command");
 
   // copy 1st arg back into line (copy is being used)
   // check maxline since arg[0] could have been exanded by variables
   // substitute for $ variables (no printing) and print arg
 
   int n = strlen(arg[0]) + 1;
   if (n > maxline) reallocate(line,maxline,n);
   strcpy(line,arg[0]);
   substitute(line,work,maxline,maxwork,0);
 
   if (me == 0) {
     if (screen) fprintf(screen,"%s\n",line);
     if (logfile) fprintf(logfile,"%s\n",line);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::quit()
 {
   if (narg) error->all(FLERR,"Illegal quit command");
   error->done();
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::shell()
 {
   if (narg < 1) error->all(FLERR,"Illegal shell command");
 
   if (strcmp(arg[0],"cd") == 0) {
     if (narg != 2) error->all(FLERR,"Illegal shell command");
     chdir(arg[1]);
 
   } else if (strcmp(arg[0],"mkdir") == 0) {
     if (narg < 2) error->all(FLERR,"Illegal shell command");
 #if !defined(WINDOWS) && !defined(__MINGW32_VERSION)
     if (me == 0)
       for (int i = 1; i < narg; i++)
         mkdir(arg[i], S_IRWXU | S_IRGRP | S_IXGRP);
 #endif
 
   } else if (strcmp(arg[0],"mv") == 0) {
     if (narg != 3) error->all(FLERR,"Illegal shell command");
     if (me == 0) rename(arg[1],arg[2]);
 
   } else if (strcmp(arg[0],"rm") == 0) {
     if (narg < 2) error->all(FLERR,"Illegal shell command");
     if (me == 0)
       for (int i = 1; i < narg; i++) unlink(arg[i]);
 
   } else if (strcmp(arg[0],"rmdir") == 0) {
     if (narg < 2) error->all(FLERR,"Illegal shell command");
     if (me == 0)
       for (int i = 1; i < narg; i++) rmdir(arg[i]);
 
   // use work string to concat args back into one string separated by spaces
   // invoke string in shell via system()
 
   } else {
     int n = 0;
     for (int i = 0; i < narg; i++) n += strlen(arg[i]) + 1;
     if (n > maxwork) reallocate(work,maxwork,n);
 
     strcpy(work,arg[0]);
     for (int i = 1; i < narg; i++) {
       strcat(work," ");
       strcat(work,arg[i]);
     }
 
     if (me == 0) system(work);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::variable_command()
 {
   variable->set(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 /* ---------------------------------------------------------------------- */
 /* ---------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
    one function for each LAMMPS-specific input script command
 ------------------------------------------------------------------------- */
 
 /* ---------------------------------------------------------------------- */
 
 void Input::angle_coeff()
 {
   if (domain->box_exist == 0)
     error->all(FLERR,"Angle_coeff command before simulation box is defined");
   if (force->angle == NULL)
     error->all(FLERR,"Angle_coeff command before angle_style is defined");
   if (atom->avec->angles_allow == 0)
     error->all(FLERR,"Angle_coeff command when no angles allowed");
   force->angle->coeff(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::angle_style()
 {
   if (narg < 1) error->all(FLERR,"Illegal angle_style command");
   if (atom->avec->angles_allow == 0)
     error->all(FLERR,"Angle_style command when no angles allowed");
   force->create_angle(arg[0],lmp->suffix);
   if (force->angle) force->angle->settings(narg-1,&arg[1]);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::atom_modify()
 {
   atom->modify_params(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::atom_style()
 {
   if (narg < 1) error->all(FLERR,"Illegal atom_style command");
   if (domain->box_exist)
     error->all(FLERR,"Atom_style command after simulation box is defined");
   atom->create_avec(arg[0],narg-1,&arg[1],lmp->suffix);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::bond_coeff()
 {
   if (domain->box_exist == 0)
     error->all(FLERR,"Bond_coeff command before simulation box is defined");
   if (force->bond == NULL)
     error->all(FLERR,"Bond_coeff command before bond_style is defined");
   if (atom->avec->bonds_allow == 0)
     error->all(FLERR,"Bond_coeff command when no bonds allowed");
   force->bond->coeff(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::bond_style()
 {
   if (narg < 1) error->all(FLERR,"Illegal bond_style command");
   if (atom->avec->bonds_allow == 0)
     error->all(FLERR,"Bond_style command when no bonds allowed");
   force->create_bond(arg[0],lmp->suffix);
   if (force->bond) force->bond->settings(narg-1,&arg[1]);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::boundary()
 {
   if (domain->box_exist)
     error->all(FLERR,"Boundary command after simulation box is defined");
   domain->set_boundary(narg,arg,0);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::box()
 {
   if (domain->box_exist)
     error->all(FLERR,"Box command after simulation box is defined");
   domain->set_box(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::communicate()
 {
   comm->set(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::compute()
 {
   modify->add_compute(narg,arg,lmp->suffix);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::compute_modify()
 {
   modify->modify_compute(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::dielectric()
 {
   if (narg != 1) error->all(FLERR,"Illegal dielectric command");
   force->dielectric = atof(arg[0]);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::dihedral_coeff()
 {
   if (domain->box_exist == 0)
     error->all(FLERR,"Dihedral_coeff command before simulation box is defined");
   if (force->dihedral == NULL)
     error->all(FLERR,"Dihedral_coeff command before dihedral_style is defined");
   if (atom->avec->dihedrals_allow == 0)
     error->all(FLERR,"Dihedral_coeff command when no dihedrals allowed");
   force->dihedral->coeff(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::dihedral_style()
 {
   if (narg < 1) error->all(FLERR,"Illegal dihedral_style command");
   if (atom->avec->dihedrals_allow == 0)
     error->all(FLERR,"Dihedral_style command when no dihedrals allowed");
   force->create_dihedral(arg[0],lmp->suffix);
   if (force->dihedral) force->dihedral->settings(narg-1,&arg[1]);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::dimension()
 {
   if (narg != 1) error->all(FLERR,"Illegal dimension command");
   if (domain->box_exist)
     error->all(FLERR,"Dimension command after simulation box is defined");
   domain->dimension = atoi(arg[0]);
   if (domain->dimension != 2 && domain->dimension != 3)
     error->all(FLERR,"Illegal dimension command");
 
   // must reset default extra_dof of all computes
   // since some were created before dimension command is encountered
 
   for (int i = 0; i < modify->ncompute; i++)
     modify->compute[i]->reset_extra_dof();
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::dump()
 {
   output->add_dump(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::dump_modify()
 {
   output->modify_dump(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::fix()
 {
   modify->add_fix(narg,arg,lmp->suffix);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::fix_modify()
 {
   modify->modify_fix(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::group_command()
 {
   group->assign(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::improper_coeff()
 {
   if (domain->box_exist == 0)
     error->all(FLERR,"Improper_coeff command before simulation box is defined");
   if (force->improper == NULL)
     error->all(FLERR,"Improper_coeff command before improper_style is defined");
   if (atom->avec->impropers_allow == 0)
     error->all(FLERR,"Improper_coeff command when no impropers allowed");
   force->improper->coeff(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::improper_style()
 {
   if (narg < 1) error->all(FLERR,"Illegal improper_style command");
   if (atom->avec->impropers_allow == 0)
     error->all(FLERR,"Improper_style command when no impropers allowed");
   force->create_improper(arg[0],lmp->suffix);
   if (force->improper) force->improper->settings(narg-1,&arg[1]);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::kspace_modify()
 {
   if (force->kspace == NULL)
     error->all(FLERR,"KSpace style has not yet been set");
   force->kspace->modify_params(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::kspace_style()
 {
   force->create_kspace(narg,arg,lmp->suffix);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::lattice()
 {
   domain->set_lattice(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::mass()
 {
   if (narg != 2) error->all(FLERR,"Illegal mass command");
   if (domain->box_exist == 0)
     error->all(FLERR,"Mass command before simulation box is defined");
   atom->set_mass(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::min_modify()
 {
   update->minimize->modify_params(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::min_style()
 {
   if (domain->box_exist == 0)
     error->all(FLERR,"Min_style command before simulation box is defined");
   update->create_minimize(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::neigh_modify()
 {
   neighbor->modify_params(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::neighbor_command()
 {
   neighbor->set(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::newton()
 {
   int newton_pair,newton_bond;
 
   if (narg == 1) {
     if (strcmp(arg[0],"off") == 0) newton_pair = newton_bond = 0;
     else if (strcmp(arg[0],"on") == 0) newton_pair = newton_bond = 1;
     else error->all(FLERR,"Illegal newton command");
   } else if (narg == 2) {
     if (strcmp(arg[0],"off") == 0) newton_pair = 0;
     else if (strcmp(arg[0],"on") == 0) newton_pair= 1;
     else error->all(FLERR,"Illegal newton command");
     if (strcmp(arg[1],"off") == 0) newton_bond = 0;
     else if (strcmp(arg[1],"on") == 0) newton_bond = 1;
     else error->all(FLERR,"Illegal newton command");
   } else error->all(FLERR,"Illegal newton command");
 
   force->newton_pair = newton_pair;
 
   if (newton_bond == 0) {
     if (domain->box_exist && force->newton_bond == 1)
       error->all(FLERR,"Newton bond change after simulation box is defined");
     force->newton_bond = 0;
   } else {
     if (domain->box_exist && force->newton_bond == 0)
       error->all(FLERR,"Newton bond change after simulation box is defined");
     force->newton_bond = 1;
   }
 
   if (newton_pair || newton_bond) force->newton = 1;
   else force->newton = 0;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::package()
 {
   if (domain->box_exist)
     error->all(FLERR,"Package command after simulation box is defined");
   if (narg < 1) error->all(FLERR,"Illegal package command");
 
   if (strcmp(arg[0],"cuda") == 0) {
     if (!lmp->cuda)
       error->all(FLERR,"Package cuda command without USER-CUDA installed");
     lmp->cuda->accelerator(narg-1,&arg[1]);
 
   } else if (strcmp(arg[0],"gpu") == 0) {
     char **fixarg = new char*[2+narg];
     fixarg[0] = (char *) "package_gpu";
     fixarg[1] = (char *) "all";
     fixarg[2] = (char *) "GPU";
     for (int i = 1; i < narg; i++) fixarg[i+2] = arg[i];
     modify->add_fix(2+narg,fixarg,NULL);
     delete [] fixarg;
     force->newton_pair = 0;
 
   } else if (strcmp(arg[0],"omp") == 0) {
     char **fixarg = new char*[2+narg];
     fixarg[0] = (char *) "package_omp";
     fixarg[1] = (char *) "all";
     fixarg[2] = (char *) "OMP";
     for (int i = 1; i < narg; i++) fixarg[i+2] = arg[i];
     modify->add_fix(2+narg,fixarg,NULL);
     delete [] fixarg;
 
   } else error->all(FLERR,"Illegal package command");
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::pair_coeff()
 {
   if (domain->box_exist == 0)
     error->all(FLERR,"Pair_coeff command before simulation box is defined");
   if (force->pair == NULL)
     error->all(FLERR,"Pair_coeff command before pair_style is defined");
   force->pair->coeff(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::pair_modify()
 {
   if (force->pair == NULL)
     error->all(FLERR,"Pair_modify command before pair_style is defined");
   force->pair->modify_params(narg,arg);
 }
 
 /* ----------------------------------------------------------------------
    if old pair style exists and new style is same, just change settings
    else create new pair class
 ------------------------------------------------------------------------- */
 
 void Input::pair_style()
 {
   if (narg < 1) error->all(FLERR,"Illegal pair_style command");
   if (force->pair && strcmp(arg[0],force->pair_style) == 0) {
     force->pair->settings(narg-1,&arg[1]);
     return;
   }
   force->create_pair(arg[0],lmp->suffix);
   if (force->pair) force->pair->settings(narg-1,&arg[1]);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::pair_write()
 {
   if (force->pair == NULL)
     error->all(FLERR,"Pair_write command before pair_style is defined");
   force->pair->write_file(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::processors()
 {
   if (domain->box_exist)
     error->all(FLERR,"Processors command after simulation box is defined");
   comm->set_processors(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::region()
 {
   domain->add_region(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::reset_timestep()
 {
   update->reset_timestep(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::restart()
 {
   output->create_restart(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::run_style()
 {
   if (domain->box_exist == 0)
     error->all(FLERR,"Run_style command before simulation box is defined");
   update->create_integrate(narg,arg,lmp->suffix);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::special_bonds()
 {
   // store 1-3,1-4 and dihedral/extra flag values before change
   // change in 1-2 coeffs will not change the special list
 
   double lj2 = force->special_lj[2];
   double lj3 = force->special_lj[3];
   double coul2 = force->special_coul[2];
   double coul3 = force->special_coul[3];
   int angle = force->special_angle;
   int dihedral = force->special_dihedral;
   int extra = force->special_extra;
 
   force->set_special(narg,arg);
 
   // if simulation box defined and saved values changed, redo special list
 
   if (domain->box_exist && atom->molecular) {
     if (lj2 != force->special_lj[2] || lj3 != force->special_lj[3] ||
         coul2 != force->special_coul[2] || coul3 != force->special_coul[3] ||
         angle != force->special_angle ||
         dihedral != force->special_dihedral ||
         extra != force->special_extra) {
       Special special(lmp);
       special.build();
     }
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::suffix()
 {
   if (narg != 1) error->all(FLERR,"Illegal suffix command");
 
   if (strcmp(arg[0],"off") == 0) lmp->suffix_enable = 0;
   else if (strcmp(arg[0],"on") == 0) lmp->suffix_enable = 1;
   else {
     delete [] lmp->suffix;
     int n = strlen(arg[0]) + 1;
     lmp->suffix = new char[n];
     strcpy(lmp->suffix,arg[0]);
     lmp->suffix_enable = 1;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::thermo()
 {
   output->set_thermo(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::thermo_modify()
 {
   output->thermo->modify_params(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::thermo_style()
 {
   output->create_thermo(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::timestep()
 {
   if (narg != 1) error->all(FLERR,"Illegal timestep command");
   update->dt = atof(arg[0]);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::uncompute()
 {
   if (narg != 1) error->all(FLERR,"Illegal uncompute command");
   modify->delete_compute(arg[0]);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::undump()
 {
   if (narg != 1) error->all(FLERR,"Illegal undump command");
   output->delete_dump(arg[0]);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::unfix()
 {
   if (narg != 1) error->all(FLERR,"Illegal unfix command");
   modify->delete_fix(arg[0]);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Input::units()
 {
   if (narg != 1) error->all(FLERR,"Illegal units command");
   if (domain->box_exist)
     error->all(FLERR,"Units command after simulation box is defined");
   update->set_units(arg[0]);
 }