diff --git a/src/Makefile.lib b/src/Makefile.lib
index 8b895ecfc..09d66d2d2 100644
--- a/src/Makefile.lib
+++ b/src/Makefile.lib
@@ -1,37 +1,63 @@
-# LAMMPS static library multiple-machine Makefile
+# LAMMPS static library multiple-machine Makefile -*- Makefile -*-
 
 SHELL = /bin/sh
 
 # Definitions
 
 ROOT =	lammps
 EXE =	lib$(ROOT)_$@.a
 
-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_template.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 citeme.cpp comm.cpp comm_brick.cpp comm_tiled.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_dilatation_atom.cpp compute_displace_atom.cpp compute_erotate_asphere.cpp compute_erotate_rigid.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_ke_rigid.cpp compute_msd.cpp compute_msd_molecule.cpp compute_msd_nongauss.cpp compute_pair.cpp compute_pair_local.cpp compute_pe.cpp compute_pe_atom.cpp compute_plasticity_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 compute_vacf.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_atom_mpiio.cpp dump_cfg.cpp dump_custom.cpp dump_custom_mpiio.cpp dump_dcd.cpp dump_image.cpp dump_local.cpp dump_movie.cpp dump_xtc.cpp dump_xyz.cpp dump_xyz_mpiio.cpp error.cpp ewald.cpp ewald_disp.cpp fft3d.cpp fft3d_wrap.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_gld.cpp fix_gravity.cpp fix_group.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_oneway.cpp fix_orient_fcc.cpp fix_peri_neigh.cpp fix_planeforce.cpp fix_pour.cpp fix_press_berendsen.cpp fix_print.cpp fix_property_atom.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_nh_small.cpp fix_rigid_nph.cpp fix_rigid_nph_small.cpp fix_rigid_npt.cpp fix_rigid_npt_small.cpp fix_rigid_nve.cpp fix_rigid_nve_small.cpp fix_rigid_nvt.cpp fix_rigid_nvt_small.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.cpp fix_store_force.cpp fix_store_state.cpp fix_temp_berendsen.cpp fix_temp_csvr.cpp fix_temp_rescale.cpp fix_thermal_conductivity.cpp fix_tmd.cpp fix_ttm.cpp fix_tune_kspace.cpp fix_vector.cpp fix_viscosity.cpp fix_viscous.cpp fix_wall.cpp fix_wall_colloid.cpp fix_wall_gran.cpp fix_wall_harmonic.cpp fix_wall_lj1043.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 gridcomm.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 molecule.cpp msm.cpp msm_cg.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_long.cpp pair_born_coul_msm.cpp pair_born_coul_wolf.cpp pair_brownian.cpp pair_brownian_poly.cpp pair_buck.cpp pair_buck_coul_cut.cpp pair_buck_coul_long.cpp pair_buck_coul_msm.cpp pair_buck_long_coul_long.cpp pair_colloid.cpp pair_comb.cpp pair_comb3.cpp pair_coul_cut.cpp pair_coul_debye.cpp pair_coul_dsf.cpp pair_coul_long.cpp pair_coul_msm.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_charmm_coul_long.cpp pair_lj_charmm_coul_long_opt.cpp pair_lj_charmm_coul_msm.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_coul_long.cpp pair_lj_cut_coul_long_opt.cpp pair_lj_cut_coul_msm.cpp pair_lj_cut_dipole_cut.cpp pair_lj_cut_dipole_long.cpp pair_lj_cut_opt.cpp pair_lj_cut_tip4p_cut.cpp pair_lj_cut_tip4p_long.cpp pair_lj_cut_tip4p_long_opt.cpp pair_lj_expand.cpp pair_lj_gromacs.cpp pair_lj_gromacs_coul_gromacs.cpp pair_lj_long_coul_long.cpp pair_lj_long_coul_long_opt.cpp pair_lj_long_dipole_long.cpp pair_lj_long_tip4p_long.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_meam.cpp pair_mie_cut.cpp pair_morse.cpp pair_morse_opt.cpp pair_nb3b_harmonic.cpp pair_nm_cut.cpp pair_nm_cut_coul_cut.cpp pair_nm_cut_coul_long.cpp pair_peri_eps.cpp pair_peri_lps.cpp pair_peri_pmb.cpp pair_peri_ves.cpp pair_rebo.cpp pair_resquared.cpp pair_soft.cpp pair_sw.cpp pair_table.cpp pair_tersoff.cpp pair_tersoff_mod.cpp pair_tersoff_zbl.cpp pair_tip4p_cut.cpp pair_tip4p_long.cpp pair_tri_lj.cpp pair_yukawa.cpp pair_yukawa_colloid.cpp pair_zbl.cpp pppm.cpp pppm_cg.cpp pppm_disp.cpp pppm_disp_tip4p.cpp pppm_stagger.cpp pppm_tip4p.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 remap.cpp remap_wrap.cpp replicate.cpp rerun.cpp respa.cpp restart_mpiio.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_data.cpp write_dump.cpp write_restart.cpp xdr_compat.cpp 
+SRC =	angle.cpp angle_charmm.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_ellipsoid.cpp atom_vec_full.cpp atom_vec_hybrid.cpp atom_vec_line.cpp atom_vec_molecular.cpp atom_vec_sphere.cpp atom_vec_template.cpp atom_vec_tri.cpp balance.cpp body.cpp bond.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 citeme.cpp comm.cpp comm_brick.cpp comm_tiled.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_dihedral_local.cpp compute_displace_atom.cpp compute_erotate_rigid.cpp compute_erotate_sphere.cpp compute_erotate_sphere_atom.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_ke_rigid.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_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_vacf.cpp create_atoms.cpp create_box.cpp delete_atoms.cpp delete_bonds.cpp dihedral.cpp dihedral_charmm.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_movie.cpp dump_xyz.cpp error.cpp finish.cpp fix.cpp fix_adapt.cpp fix_addforce.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_box_relax.cpp fix_deform.cpp fix_drag.cpp fix_dt_reset.cpp fix_enforce2d.cpp fix_external.cpp fix_freeze.cpp fix_gravity.cpp fix_group.cpp fix_heat.cpp fix_indent.cpp fix_langevin.cpp fix_lineforce.cpp fix_minimize.cpp fix_momentum.cpp fix_move.cpp fix_nh.cpp fix_nh_sphere.cpp fix_nph.cpp fix_nph_sphere.cpp fix_npt.cpp fix_npt_sphere.cpp fix_nve.cpp fix_nve_limit.cpp fix_nve_noforce.cpp fix_nve_sphere.cpp fix_nvt.cpp fix_nvt_sllod.cpp fix_nvt_sphere.cpp fix_planeforce.cpp fix_pour.cpp fix_press_berendsen.cpp fix_print.cpp fix_property_atom.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_nh_small.cpp fix_rigid_nph.cpp fix_rigid_nph_small.cpp fix_rigid_npt.cpp fix_rigid_npt_small.cpp fix_rigid_nve.cpp fix_rigid_nve_small.cpp fix_rigid_nvt.cpp fix_rigid_nvt_small.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_store.cpp fix_store_force.cpp fix_store_state.cpp fix_temp_berendsen.cpp fix_temp_csvr.cpp fix_temp_rescale.cpp fix_tmd.cpp fix_vector.cpp fix_viscous.cpp fix_wall.cpp fix_wall_gran.cpp fix_wall_harmonic.cpp fix_wall_lj1043.cpp fix_wall_lj126.cpp fix_wall_lj93.cpp fix_wall_reflect.cpp fix_wall_region.cpp force.cpp group.cpp image.cpp improper.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 molecule.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_buck.cpp pair_buck_coul_cut.cpp pair_comb.cpp pair_comb3.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_eam.cpp pair_eam_alloy.cpp pair_eam_fs.cpp pair_eim.cpp pair_gauss.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_lj96_cut.cpp pair_lj_charmm_coul_charmm.cpp pair_lj_charmm_coul_charmm_implicit.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_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_mie_cut.cpp pair_morse.cpp pair_nb3b_harmonic.cpp pair_rebo.cpp pair_soft.cpp pair_sw.cpp pair_table.cpp pair_tersoff.cpp pair_tersoff_mod.cpp pair_tersoff_zbl.cpp pair_tip4p_cut.cpp pair_yukawa.cpp pair_zbl.cpp procmap.cpp random_mars.cpp random_park.cpp rcb.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 thermo.cpp timer.cpp universe.cpp update.cpp variable.cpp velocity.cpp verlet.cpp write_data.cpp write_dump.cpp write_restart.cpp 
 
-INC =	accelerator_cuda.h accelerator_kokkos.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_template.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 citeme.h comm.h comm_brick.h comm_tiled.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_dilatation_atom.h compute_displace_atom.h compute_erotate_asphere.h compute_erotate_rigid.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_ke_rigid.h compute_msd.h compute_msd_molecule.h compute_msd_nongauss.h compute_pair.h compute_pair_local.h compute_pe.h compute_pe_atom.h compute_plasticity_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 compute_vacf.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_atom_mpiio.h dump_cfg.h dump_custom.h dump_custom_mpiio.h dump_dcd.h dump_image.h dump_local.h dump_movie.h dump_xtc.h dump_xyz.h dump_xyz_mpiio.h error.h ewald.h ewald_disp.h fft3d.h fft3d_wrap.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_gld.h fix_gravity.h fix_group.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_oneway.h fix_orient_fcc.h fix_peri_neigh.h fix_planeforce.h fix_pour.h fix_press_berendsen.h fix_print.h fix_property_atom.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_nh_small.h fix_rigid_nph.h fix_rigid_nph_small.h fix_rigid_npt.h fix_rigid_npt_small.h fix_rigid_nve.h fix_rigid_nve_small.h fix_rigid_nvt.h fix_rigid_nvt_small.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.h fix_store_force.h fix_store_state.h fix_temp_berendsen.h fix_temp_csvr.h fix_temp_rescale.h fix_thermal_conductivity.h fix_tmd.h fix_ttm.h fix_tune_kspace.h fix_vector.h fix_viscosity.h fix_viscous.h fix_wall.h fix_wall_colloid.h fix_wall_gran.h fix_wall_harmonic.h fix_wall_lj1043.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 gridcomm.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 kissfft.h kspace.h lammps.h lattice.h library.h lmptype.h lmpwindows.h math_complex.h math_const.h math_extra.h math_special.h math_vector.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 molecule.h mpiio.h msm.h msm_cg.h my_page.h my_pool_chunk.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_long.h pair_born_coul_msm.h pair_born_coul_wolf.h pair_brownian.h pair_brownian_poly.h pair_buck.h pair_buck_coul_cut.h pair_buck_coul_long.h pair_buck_coul_msm.h pair_buck_long_coul_long.h pair_colloid.h pair_comb.h pair_comb3.h pair_coul_cut.h pair_coul_debye.h pair_coul_dsf.h pair_coul_long.h pair_coul_msm.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_charmm_coul_long.h pair_lj_charmm_coul_long_opt.h pair_lj_charmm_coul_msm.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_coul_long.h pair_lj_cut_coul_long_opt.h pair_lj_cut_coul_msm.h pair_lj_cut_dipole_cut.h pair_lj_cut_dipole_long.h pair_lj_cut_opt.h pair_lj_cut_tip4p_cut.h pair_lj_cut_tip4p_long.h pair_lj_cut_tip4p_long_opt.h pair_lj_expand.h pair_lj_gromacs.h pair_lj_gromacs_coul_gromacs.h pair_lj_long_coul_long.h pair_lj_long_coul_long_opt.h pair_lj_long_dipole_long.h pair_lj_long_tip4p_long.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_meam.h pair_mie_cut.h pair_morse.h pair_morse_opt.h pair_nb3b_harmonic.h pair_nm_cut.h pair_nm_cut_coul_cut.h pair_nm_cut_coul_long.h pair_peri_eps.h pair_peri_lps.h pair_peri_pmb.h pair_peri_ves.h pair_rebo.h pair_resquared.h pair_soft.h pair_sw.h pair_table.h pair_tersoff.h pair_tersoff_mod.h pair_tersoff_zbl.h pair_tip4p_cut.h pair_tip4p_long.h pair_tri_lj.h pair_yukawa.h pair_yukawa_colloid.h pair_zbl.h pointers.h pppm.h pppm_cg.h pppm_disp.h pppm_disp_tip4p.h pppm_stagger.h pppm_tip4p.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 remap.h remap_wrap.h replicate.h rerun.h respa.h restart_mpiio.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_data.h write_dump.h write_restart.h xdr_compat.h 
+INC =	accelerator_cuda.h accelerator_intel.h accelerator_kokkos.h accelerator_omp.h angle.h angle_charmm.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_ellipsoid.h atom_vec_full.h atom_vec_hybrid.h atom_vec_line.h atom_vec_molecular.h atom_vec_sphere.h atom_vec_template.h atom_vec_tri.h balance.h body.h bond.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 citeme.h comm.h comm_brick.h comm_tiled.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_dihedral_local.h compute_displace_atom.h compute_erotate_rigid.h compute_erotate_sphere.h compute_erotate_sphere_atom.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_ke_rigid.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_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_vacf.h create_atoms.h create_box.h delete_atoms.h delete_bonds.h dihedral.h dihedral_charmm.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_movie.h dump_xyz.h error.h finish.h fix.h fix_adapt.h fix_addforce.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_box_relax.h fix_deform.h fix_drag.h fix_dt_reset.h fix_enforce2d.h fix_external.h fix_freeze.h fix_gravity.h fix_group.h fix_heat.h fix_indent.h fix_langevin.h fix_lineforce.h fix_minimize.h fix_momentum.h fix_move.h fix_nh.h fix_nh_sphere.h fix_nph.h fix_nph_sphere.h fix_npt.h fix_npt_sphere.h fix_nve.h fix_nve_limit.h fix_nve_noforce.h fix_nve_sphere.h fix_nvt.h fix_nvt_sllod.h fix_nvt_sphere.h fix_planeforce.h fix_pour.h fix_press_berendsen.h fix_print.h fix_property_atom.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_nh_small.h fix_rigid_nph.h fix_rigid_nph_small.h fix_rigid_npt.h fix_rigid_npt_small.h fix_rigid_nve.h fix_rigid_nve_small.h fix_rigid_nvt.h fix_rigid_nvt_small.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_store.h fix_store_force.h fix_store_state.h fix_temp_berendsen.h fix_temp_csvr.h fix_temp_rescale.h fix_tmd.h fix_vector.h fix_viscous.h fix_wall.h fix_wall_gran.h fix_wall_harmonic.h fix_wall_lj1043.h fix_wall_lj126.h fix_wall_lj93.h fix_wall_reflect.h fix_wall_region.h force.h group.h image.h improper.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_complex.h math_const.h math_extra.h math_special.h math_vector.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 molecule.h mpiio.h my_page.h my_pool_chunk.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_buck.h pair_buck_coul_cut.h pair_comb.h pair_comb3.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_eam.h pair_eam_alloy.h pair_eam_fs.h pair_eim.h pair_gauss.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_lj96_cut.h pair_lj_charmm_coul_charmm.h pair_lj_charmm_coul_charmm_implicit.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_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_mie_cut.h pair_morse.h pair_nb3b_harmonic.h pair_rebo.h pair_soft.h pair_sw.h pair_table.h pair_tersoff.h pair_tersoff_mod.h pair_tersoff_zbl.h pair_tip4p_cut.h pair_yukawa.h pair_zbl.h pointers.h procmap.h random_mars.h random_park.h rcb.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 thermo.h timer.h universe.h update.h variable.h velocity.h verlet.h version.h write_data.h write_dump.h write_restart.h 
 
 OBJ =	$(SRC:.cpp=.o)
 
 # Targets
 
 help:
-	@echo 'Type "make target" where target is one of:'
+	@echo 'Type "make -f Makefile.lib target" where target is one of:'
 	@echo ''
 	@files="`ls MAKE/Makefile.*`"; \
-	for file in $$files; do head -1 $$file; done
+	  for file in $$files; do head -1 $$file; done
+	@echo ''
+	@echo 'or one of these from src/MAKE/OPTIONS:'
+	@echo ''
+	@files="`ls MAKE/OPTIONS/Makefile.*`"; \
+	  for file in $$files; do head -1 $$file; done
+	@echo ''
+	@echo 'or one of these from src/MAKE/MACHINES:'
+	@echo ''
+	@files="`ls MAKE/MACHINES/Makefile.*`"; \
+	  for file in $$files; do head -1 $$file; done
+	@echo ''
+	@echo 'or one of these from src/MAKE/MINE:'
+	@echo ''
+	@files="`ls MAKE/MINE/Makefile.* 2>/dev/null`"; \
+	  for file in $$files; do head -1 $$file; done
+	@echo ''
 
 .DEFAULT:
-	@test -f MAKE/Makefile.$@
+	@if [ $@ = "serial" -a ! -f STUBS/libmpi_stubs.a ]; \
+	  then $(MAKE) stubs; fi
+	@test -f MAKE/Makefile.$@ -o -f MAKE/OPTIONS/Makefile.$@ -o \
+	  -f MAKE/MACHINES/Makefile.$@ -o -f MAKE/MINE/Makefile.$@
 	@if [ ! -d Obj_$@ ]; then mkdir Obj_$@; fi
-	@cp -p $(SRC) $(INC) Obj_$@
-	@cp MAKE/Makefile.$@ Obj_$@/Makefile
+	@$(SHELL) Make.sh style
+	@if [ -f MAKE/MACHINES/Makefile.$@ ]; \
+	  then cp MAKE/MACHINES/Makefile.$@ Obj_$@/Makefile; fi
+	@if [ -f MAKE/OPTIONS/Makefile.$@ ]; \
+	  then cp MAKE/OPTIONS/Makefile.$@ Obj_$@/Makefile; fi
+	@if [ -f MAKE/Makefile.$@ ]; \
+	  then cp MAKE/Makefile.$@ Obj_$@/Makefile; fi
+	@if [ -f MAKE/MINE/Makefile.$@ ]; \
+	  then cp MAKE/MINE/Makefile.$@ Obj_$@/Makefile; fi
 	@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_$@
 	@cd Obj_$@; \
 	$(MAKE) $(MFLAGS) "OBJ = $(OBJ)" "INC = $(INC)" "SHFLAGS =" \
           "EXE = ../$(EXE)" lib
 	@if [ -d Obj_$@ ]; then cd Obj_$@; rm -f $(SRC) $(INC) Makefile*; fi
diff --git a/src/Makefile.shlib b/src/Makefile.shlib
index fd00ed128..d14216466 100644
--- a/src/Makefile.shlib
+++ b/src/Makefile.shlib
@@ -1,80 +1,80 @@
 # 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_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_template.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 citeme.cpp comm.cpp comm_brick.cpp comm_tiled.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_dilatation_atom.cpp compute_displace_atom.cpp compute_erotate_asphere.cpp compute_erotate_rigid.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_ke_rigid.cpp compute_msd.cpp compute_msd_molecule.cpp compute_msd_nongauss.cpp compute_pair.cpp compute_pair_local.cpp compute_pe.cpp compute_pe_atom.cpp compute_plasticity_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 compute_vacf.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_movie.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_gld.cpp fix_gravity.cpp fix_group.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_oneway.cpp fix_orient_fcc.cpp fix_peri_neigh.cpp fix_planeforce.cpp fix_pour.cpp fix_press_berendsen.cpp fix_print.cpp fix_property_atom.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_nh_small.cpp fix_rigid_nph.cpp fix_rigid_nph_small.cpp fix_rigid_npt.cpp fix_rigid_npt_small.cpp fix_rigid_nve.cpp fix_rigid_nve_small.cpp fix_rigid_nvt.cpp fix_rigid_nvt_small.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.cpp fix_store_force.cpp fix_store_state.cpp fix_temp_berendsen.cpp fix_temp_csvr.cpp fix_temp_rescale.cpp fix_thermal_conductivity.cpp fix_tmd.cpp fix_ttm.cpp fix_vector.cpp fix_viscosity.cpp fix_viscous.cpp fix_wall.cpp fix_wall_colloid.cpp fix_wall_gran.cpp fix_wall_harmonic.cpp fix_wall_lj1043.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 molecule.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_comb3.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_dipole_cut.cpp pair_lj_cut_dipole_long.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_long_dipole_long.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_nb3b_harmonic.cpp pair_nm_cut.cpp pair_nm_cut_coul_cut.cpp pair_nm_cut_coul_long.cpp pair_peri_eps.cpp pair_peri_lps.cpp pair_peri_pmb.cpp pair_peri_ves.cpp pair_rebo.cpp pair_resquared.cpp pair_soft.cpp pair_sw.cpp pair_table.cpp pair_tersoff.cpp pair_tersoff_mod.cpp pair_tersoff_zbl.cpp pair_tip4p_cut.cpp pair_tri_lj.cpp pair_yukawa.cpp pair_yukawa_colloid.cpp pair_zbl.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_data.cpp write_dump.cpp write_restart.cpp xdr_compat.cpp 
+SRC =	angle.cpp angle_charmm.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_ellipsoid.cpp atom_vec_full.cpp atom_vec_hybrid.cpp atom_vec_line.cpp atom_vec_molecular.cpp atom_vec_sphere.cpp atom_vec_template.cpp atom_vec_tri.cpp balance.cpp body.cpp bond.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 citeme.cpp comm.cpp comm_brick.cpp comm_tiled.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_dihedral_local.cpp compute_displace_atom.cpp compute_erotate_rigid.cpp compute_erotate_sphere.cpp compute_erotate_sphere_atom.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_ke_rigid.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_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_vacf.cpp create_atoms.cpp create_box.cpp delete_atoms.cpp delete_bonds.cpp dihedral.cpp dihedral_charmm.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_movie.cpp dump_xyz.cpp error.cpp finish.cpp fix.cpp fix_adapt.cpp fix_addforce.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_box_relax.cpp fix_deform.cpp fix_drag.cpp fix_dt_reset.cpp fix_enforce2d.cpp fix_external.cpp fix_freeze.cpp fix_gravity.cpp fix_group.cpp fix_heat.cpp fix_indent.cpp fix_langevin.cpp fix_lineforce.cpp fix_minimize.cpp fix_momentum.cpp fix_move.cpp fix_nh.cpp fix_nh_sphere.cpp fix_nph.cpp fix_nph_sphere.cpp fix_npt.cpp fix_npt_sphere.cpp fix_nve.cpp fix_nve_limit.cpp fix_nve_noforce.cpp fix_nve_sphere.cpp fix_nvt.cpp fix_nvt_sllod.cpp fix_nvt_sphere.cpp fix_planeforce.cpp fix_pour.cpp fix_press_berendsen.cpp fix_print.cpp fix_property_atom.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_nh_small.cpp fix_rigid_nph.cpp fix_rigid_nph_small.cpp fix_rigid_npt.cpp fix_rigid_npt_small.cpp fix_rigid_nve.cpp fix_rigid_nve_small.cpp fix_rigid_nvt.cpp fix_rigid_nvt_small.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_store.cpp fix_store_force.cpp fix_store_state.cpp fix_temp_berendsen.cpp fix_temp_csvr.cpp fix_temp_rescale.cpp fix_tmd.cpp fix_vector.cpp fix_viscous.cpp fix_wall.cpp fix_wall_gran.cpp fix_wall_harmonic.cpp fix_wall_lj1043.cpp fix_wall_lj126.cpp fix_wall_lj93.cpp fix_wall_reflect.cpp fix_wall_region.cpp force.cpp group.cpp image.cpp improper.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 molecule.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_buck.cpp pair_buck_coul_cut.cpp pair_comb.cpp pair_comb3.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_eam.cpp pair_eam_alloy.cpp pair_eam_fs.cpp pair_eim.cpp pair_gauss.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_lj96_cut.cpp pair_lj_charmm_coul_charmm.cpp pair_lj_charmm_coul_charmm_implicit.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_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_mie_cut.cpp pair_morse.cpp pair_nb3b_harmonic.cpp pair_rebo.cpp pair_soft.cpp pair_sw.cpp pair_table.cpp pair_tersoff.cpp pair_tersoff_mod.cpp pair_tersoff_zbl.cpp pair_tip4p_cut.cpp pair_yukawa.cpp pair_zbl.cpp procmap.cpp random_mars.cpp random_park.cpp rcb.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 thermo.cpp timer.cpp universe.cpp update.cpp variable.cpp velocity.cpp verlet.cpp write_data.cpp write_dump.cpp write_restart.cpp 
 
-INC =	accelerator_cuda.h accelerator_kokkos.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_template.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 citeme.h comm.h comm_brick.h comm_tiled.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_dilatation_atom.h compute_displace_atom.h compute_erotate_asphere.h compute_erotate_rigid.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_ke_rigid.h compute_msd.h compute_msd_molecule.h compute_msd_nongauss.h compute_pair.h compute_pair_local.h compute_pe.h compute_pe_atom.h compute_plasticity_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 compute_vacf.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_movie.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_gld.h fix_gravity.h fix_group.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_oneway.h fix_orient_fcc.h fix_peri_neigh.h fix_planeforce.h fix_pour.h fix_press_berendsen.h fix_print.h fix_property_atom.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_nh_small.h fix_rigid_nph.h fix_rigid_nph_small.h fix_rigid_npt.h fix_rigid_npt_small.h fix_rigid_nve.h fix_rigid_nve_small.h fix_rigid_nvt.h fix_rigid_nvt_small.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.h fix_store_force.h fix_store_state.h fix_temp_berendsen.h fix_temp_csvr.h fix_temp_rescale.h fix_thermal_conductivity.h fix_tmd.h fix_ttm.h fix_vector.h fix_viscosity.h fix_viscous.h fix_wall.h fix_wall_colloid.h fix_wall_gran.h fix_wall_harmonic.h fix_wall_lj1043.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_complex.h math_const.h math_extra.h math_special.h math_vector.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 molecule.h mpiio.h my_page.h my_pool_chunk.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_comb3.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_dipole_cut.h pair_lj_cut_dipole_long.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_long_dipole_long.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_nb3b_harmonic.h pair_nm_cut.h pair_nm_cut_coul_cut.h pair_nm_cut_coul_long.h pair_peri_eps.h pair_peri_lps.h pair_peri_pmb.h pair_peri_ves.h pair_rebo.h pair_resquared.h pair_soft.h pair_sw.h pair_table.h pair_tersoff.h pair_tersoff_mod.h pair_tersoff_zbl.h pair_tip4p_cut.h pair_tri_lj.h pair_yukawa.h pair_yukawa_colloid.h pair_zbl.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_data.h write_dump.h write_restart.h xdr_compat.h 
+INC =	accelerator_cuda.h accelerator_intel.h accelerator_kokkos.h accelerator_omp.h angle.h angle_charmm.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_ellipsoid.h atom_vec_full.h atom_vec_hybrid.h atom_vec_line.h atom_vec_molecular.h atom_vec_sphere.h atom_vec_template.h atom_vec_tri.h balance.h body.h bond.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 citeme.h comm.h comm_brick.h comm_tiled.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_dihedral_local.h compute_displace_atom.h compute_erotate_rigid.h compute_erotate_sphere.h compute_erotate_sphere_atom.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_ke_rigid.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_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_vacf.h create_atoms.h create_box.h delete_atoms.h delete_bonds.h dihedral.h dihedral_charmm.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_movie.h dump_xyz.h error.h finish.h fix.h fix_adapt.h fix_addforce.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_box_relax.h fix_deform.h fix_drag.h fix_dt_reset.h fix_enforce2d.h fix_external.h fix_freeze.h fix_gravity.h fix_group.h fix_heat.h fix_indent.h fix_langevin.h fix_lineforce.h fix_minimize.h fix_momentum.h fix_move.h fix_nh.h fix_nh_sphere.h fix_nph.h fix_nph_sphere.h fix_npt.h fix_npt_sphere.h fix_nve.h fix_nve_limit.h fix_nve_noforce.h fix_nve_sphere.h fix_nvt.h fix_nvt_sllod.h fix_nvt_sphere.h fix_planeforce.h fix_pour.h fix_press_berendsen.h fix_print.h fix_property_atom.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_nh_small.h fix_rigid_nph.h fix_rigid_nph_small.h fix_rigid_npt.h fix_rigid_npt_small.h fix_rigid_nve.h fix_rigid_nve_small.h fix_rigid_nvt.h fix_rigid_nvt_small.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_store.h fix_store_force.h fix_store_state.h fix_temp_berendsen.h fix_temp_csvr.h fix_temp_rescale.h fix_tmd.h fix_vector.h fix_viscous.h fix_wall.h fix_wall_gran.h fix_wall_harmonic.h fix_wall_lj1043.h fix_wall_lj126.h fix_wall_lj93.h fix_wall_reflect.h fix_wall_region.h force.h group.h image.h improper.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_complex.h math_const.h math_extra.h math_special.h math_vector.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 molecule.h mpiio.h my_page.h my_pool_chunk.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_buck.h pair_buck_coul_cut.h pair_comb.h pair_comb3.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_eam.h pair_eam_alloy.h pair_eam_fs.h pair_eim.h pair_gauss.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_lj96_cut.h pair_lj_charmm_coul_charmm.h pair_lj_charmm_coul_charmm_implicit.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_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_mie_cut.h pair_morse.h pair_nb3b_harmonic.h pair_rebo.h pair_soft.h pair_sw.h pair_table.h pair_tersoff.h pair_tersoff_mod.h pair_tersoff_zbl.h pair_tip4p_cut.h pair_yukawa.h pair_zbl.h pointers.h procmap.h random_mars.h random_park.h rcb.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 thermo.h timer.h universe.h update.h variable.h velocity.h verlet.h version.h write_data.h write_dump.h write_restart.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
 	@echo ''
 	@echo 'or one of these from src/MAKE/OPTIONS:'
 	@echo ''
 	@files="`ls MAKE/OPTIONS/Makefile.*`"; \
 	  for file in $$files; do head -1 $$file; done
 	@echo ''
 	@echo 'or one of these from src/MAKE/MACHINES:'
 	@echo ''
 	@files="`ls MAKE/MACHINES/Makefile.*`"; \
 	  for file in $$files; do head -1 $$file; done
 	@echo ''
 	@echo 'or one of these from src/MAKE/MINE:'
 	@echo ''
 	@files="`ls MAKE/MINE/Makefile.* 2>/dev/null`"; \
 	  for file in $$files; do head -1 $$file; done
 	@echo ''
 
 
 .DEFAULT:
 	@if [ $@ = "serial" -a ! -f STUBS/libmpi_stubs.a ]; \
 	  then $(MAKE) stubs; fi
 	@test -f MAKE/Makefile.$@ -o -f MAKE/OPTIONS/Makefile.$@ -o \
 	  -f MAKE/MACHINES/Makefile.$@ -o -f MAKE/MINE/Makefile.$@
 	@if [ ! -d Obj_shlib_$@ ]; then mkdir Obj_shlib_$@; fi
 	@$(SHELL) Make.sh style
 	@if [ -f MAKE/MACHINES/Makefile.$@ ]; \
 	  then cp MAKE/MACHINES/Makefile.$@ Obj_shlib_$@/Makefile; fi
 	@if [ -f MAKE/OPTIONS/Makefile.$@ ]; \
 	  then cp MAKE/OPTIONS/Makefile.$@ Obj_shlib_$@/Makefile; fi
 	@if [ -f MAKE/Makefile.$@ ]; \
 	  then cp MAKE/Makefile.$@ Obj_shlib_$@/Makefile; fi
 	@if [ -f MAKE/MINE/Makefile.$@ ]; \
 	  then cp MAKE/MINE/Makefile.$@ Obj_shlib_$@/Makefile; fi
 	@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
 
 # Remove machine-specific object files
 
 clean:
 	@echo 'make clean-all           delete all object files'
 	@echo 'make clean-machine       delete object files for one machine'
 
 clean-all:
 	rm -rf Obj_shlib_*
 
 clean-%:
 	rm -rf Obj_shlib_$(@:clean-%=%)
 
diff --git a/src/compute_pair_local.cpp b/src/compute_pair_local.cpp
index 51985846a..2c19e30af 100644
--- a/src/compute_pair_local.cpp
+++ b/src/compute_pair_local.cpp
@@ -1,289 +1,289 @@
 /* ----------------------------------------------------------------------
    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 "math.h"
 #include "string.h"
 #include "stdlib.h"
 #include "compute_pair_local.h"
 #include "atom.h"
 #include "update.h"
 #include "force.h"
 #include "pair.h"
 #include "neighbor.h"
 #include "neigh_request.h"
 #include "neigh_list.h"
 #include "group.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
 #define DELTA 10000
 
 enum{DIST,ENG,FORCE,FX,FY,FZ,PN};
 
 /* ---------------------------------------------------------------------- */
 
 ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg)
 {
   if (narg < 4) error->all(FLERR,"Illegal compute pair/local command");
 
   local_flag = 1;
   nvalues = narg - 3;
   if (nvalues == 1) size_local_cols = 0;
   else size_local_cols = nvalues;
 
   pstyle = new int[nvalues];
   pindex = new int[nvalues];
 
   nvalues = 0;
   for (int iarg = 3; iarg < narg; iarg++) {
     if (strcmp(arg[iarg],"dist") == 0) pstyle[nvalues++] = DIST;
     else if (strcmp(arg[iarg],"eng") == 0) pstyle[nvalues++] = ENG;
     else if (strcmp(arg[iarg],"force") == 0) pstyle[nvalues++] = FORCE;
     else if (strcmp(arg[iarg],"fx") == 0) pstyle[nvalues++] = FX;
     else if (strcmp(arg[iarg],"fy") == 0) pstyle[nvalues++] = FY;
     else if (strcmp(arg[iarg],"fz") == 0) pstyle[nvalues++] = FZ;
     else if (arg[iarg][0] == 'p') {
       int n = atoi(&arg[iarg][1]);
       if (n <= 0) error->all(FLERR,
                              "Invalid keyword in compute pair/local command");
       pstyle[nvalues] = PN;
       pindex[nvalues++] = n-1;
     } else error->all(FLERR,"Invalid keyword in compute pair/local command");
   }
 
   // set singleflag if need to call pair->single()
 
   singleflag = 0;
   for (int i = 0; i < nvalues; i++)
     if (pstyle[i] != DIST) singleflag = 1;
 
   nmax = 0;
   vector = NULL;
   array = NULL;
 }
 
 /* ---------------------------------------------------------------------- */
 
 ComputePairLocal::~ComputePairLocal()
 {
   memory->destroy(vector);
   memory->destroy(array);
   delete [] pstyle;
   delete [] pindex;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePairLocal::init()
 {
   if (singleflag && force->pair == NULL)
     error->all(FLERR,"No pair style is defined for compute pair/local");
   if (singleflag && force->pair->single_enable == 0)
     error->all(FLERR,"Pair style does not support compute pair/local");
 
   for (int i = 0; i < nvalues; i++)
     if (pstyle[i] == PN && pindex[i] >= force->pair->single_extra)
       error->all(FLERR,"Pair style does not have extra field"
                  " requested by compute pair/local");
 
   // need an occasional half neighbor list
 
   int irequest = neighbor->request((void *) this);
   neighbor->requests[irequest]->pair = 0;
   neighbor->requests[irequest]->compute = 1;
   neighbor->requests[irequest]->occasional = 1;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePairLocal::init_list(int id, NeighList *ptr)
 {
   list = ptr;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePairLocal::compute_local()
 {
   invoked_local = update->ntimestep;
 
   // count local entries and compute pair info
 
   ncount = compute_pairs(0);
   if (ncount > nmax) reallocate(ncount);
   size_local_rows = ncount;
   compute_pairs(1);
 }
 
 /* ----------------------------------------------------------------------
    count pairs and compute pair info on this proc
    only count pair once if newton_pair is off
    both atom I,J must be in group
    if flag is set, compute requested info about pair
 ------------------------------------------------------------------------- */
 
 int ComputePairLocal::compute_pairs(int flag)
 {
   int i,j,m,n,ii,jj,inum,jnum,itype,jtype;
   tagint itag,jtag;
   double xtmp,ytmp,ztmp,delx,dely,delz;
   double rsq,eng,fpair,factor_coul,factor_lj;
   int *ilist,*jlist,*numneigh,**firstneigh;
   double *ptr;
 
   double **x = atom->x;
   tagint *tag = atom->tag;
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
   double *special_coul = force->special_coul;
   double *special_lj = force->special_lj;
   int newton_pair = force->newton_pair;
 
   // invoke half neighbor list (will copy or build if necessary)
 
   if (flag == 0) neighbor->build_one(list);
 
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
   // loop over neighbors of my atoms
   // skip if I or J are not in group
   // for newton = 0 and J = ghost atom,
   //   need to insure I,J pair is only output by one proc
-  //   use same itag,jtag lobic as in Neighbor::neigh_half_nsq()
+  //   use same itag,jtag logic as in Neighbor::neigh_half_nsq()
   // for flag = 0, just count pair interactions within force cutoff
   // for flag = 1, calculate requested output fields
 
   Pair *pair = force->pair;
   double **cutsq = force->pair->cutsq;
 
   m = 0;
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
     if (!(mask[i] & groupbit)) continue;
 
     xtmp = x[i][0];
     ytmp = x[i][1];
     ztmp = x[i][2];
     itag = tag[i];
     itype = type[i];
     jlist = firstneigh[i];
     jnum = numneigh[i];
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
       factor_lj = special_lj[sbmask(j)];
       factor_coul = special_coul[sbmask(j)];
       j &= NEIGHMASK;
 
       if (!(mask[j] & groupbit)) continue;
 
       // itag = jtag is possible for long cutoffs that include images of self
 
       if (newton_pair == 0 && j >= nlocal) {
         jtag = tag[j];
         if (itag > jtag) {
           if ((itag+jtag) % 2 == 0) continue;
         } else if (itag < jtag) {
           if ((itag+jtag) % 2 == 1) continue;
         } else {
           if (x[j][2] < ztmp) continue;
           if (x[j][2] == ztmp) {
             if (x[j][1] < ytmp) continue;
             if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
           }
         }
       }
 
       delx = xtmp - x[j][0];
       dely = ytmp - x[j][1];
       delz = ztmp - x[j][2];
       rsq = delx*delx + dely*dely + delz*delz;
       jtype = type[j];
       if (rsq >= cutsq[itype][jtype]) continue;
 
       if (flag) {
         if (singleflag)
           eng = pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair);
 
         if (nvalues == 1) ptr = &vector[m];
         else ptr = array[m];
 
         for (n = 0; n < nvalues; n++) {
           switch (pstyle[n]) {
           case DIST:
             ptr[n] = sqrt(rsq);
             break;
           case ENG:
             ptr[n] = eng;
             break;
           case FORCE:
             ptr[n] = sqrt(rsq)*fpair;
             break;
           case FX:
             ptr[n] = delx*fpair;
             break;
           case FY:
             ptr[n] = dely*fpair;
             break;
           case FZ:
             ptr[n] = delz*fpair;
             break;
           case PN:
             ptr[n] = pair->svector[pindex[n]];
             break;
           }
         }
       }
 
       m++;
     }
   }
 
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePairLocal::reallocate(int n)
 {
   // grow vector or array and indices array
 
   while (nmax < n) nmax += DELTA;
 
   if (nvalues == 1) {
     memory->destroy(vector);
     memory->create(vector,nmax,"pair/local:vector");
     vector_local = vector;
   } else {
     memory->destroy(array);
     memory->create(array,nmax,nvalues,"pair/local:array");
     array_local = array;
   }
 }
 
 /* ----------------------------------------------------------------------
    memory usage of local data
 ------------------------------------------------------------------------- */
 
 double ComputePairLocal::memory_usage()
 {
   double bytes = nmax*nvalues * sizeof(double);
   return bytes;
 }
diff --git a/src/compute_property_local.cpp b/src/compute_property_local.cpp
index d5cc0d73d..fdfd052b6 100644
--- a/src/compute_property_local.cpp
+++ b/src/compute_property_local.cpp
@@ -1,915 +1,937 @@
 /* ----------------------------------------------------------------------
    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 "string.h"
 #include "compute_property_local.h"
 #include "atom.h"
 #include "atom_vec.h"
 #include "update.h"
 #include "force.h"
 #include "pair.h"
 #include "neighbor.h"
 #include "neigh_request.h"
 #include "neigh_list.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
 enum{NONE,NEIGH,PAIR,BOND,ANGLE,DIHEDRAL,IMPROPER};
 
 #define DELTA 10000
 
 /* ---------------------------------------------------------------------- */
 
 ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg)
 {
   if (narg < 4) error->all(FLERR,"Illegal compute property/local command");
 
   local_flag = 1;
   nvalues = narg - 3;
   if (nvalues == 1) size_local_cols = 0;
   else size_local_cols = nvalues;
 
   pack_choice = new FnPtrPack[nvalues];
 
   kindflag = NONE;
 
   int i;
   for (int iarg = 3; iarg < narg; iarg++) {
     i = iarg-3;
 
     if (strcmp(arg[iarg],"natom1") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_patom1;
       if (kindflag != NONE && kindflag != NEIGH)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = NEIGH;
     } else if (strcmp(arg[iarg],"natom2") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_patom2;
       if (kindflag != NONE && kindflag != NEIGH)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = NEIGH;
     } else if (strcmp(arg[iarg],"ntype1") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_ptype1;
       if (kindflag != NONE && kindflag != NEIGH)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = NEIGH;
     } else if (strcmp(arg[iarg],"ntype2") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_ptype2;
       if (kindflag != NONE && kindflag != NEIGH)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = NEIGH;
 
     } else if (strcmp(arg[iarg],"patom1") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_patom1;
       if (kindflag != NONE && kindflag != PAIR)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = PAIR;
     } else if (strcmp(arg[iarg],"patom2") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_patom2;
       if (kindflag != NONE && kindflag != PAIR)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = PAIR;
     } else if (strcmp(arg[iarg],"ptype1") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_ptype1;
       if (kindflag != NONE && kindflag != PAIR)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = PAIR;
     } else if (strcmp(arg[iarg],"ptype2") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_ptype2;
       if (kindflag != NONE && kindflag != PAIR)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = PAIR;
 
     } else if (strcmp(arg[iarg],"batom1") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_batom1;
       if (kindflag != NONE && kindflag != BOND)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = BOND;
     } else if (strcmp(arg[iarg],"batom2") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_batom2;
       if (kindflag != NONE && kindflag != BOND)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = BOND;
     } else if (strcmp(arg[iarg],"btype") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_btype;
       if (kindflag != NONE && kindflag != BOND)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = BOND;
 
     } else if (strcmp(arg[iarg],"aatom1") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_aatom1;
       if (kindflag != NONE && kindflag != ANGLE)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = ANGLE;
     } else if (strcmp(arg[iarg],"aatom2") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_aatom2;
       if (kindflag != NONE && kindflag != ANGLE)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = ANGLE;
     } else if (strcmp(arg[iarg],"aatom3") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_aatom3;
       if (kindflag != NONE && kindflag != ANGLE)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = ANGLE;
     } else if (strcmp(arg[iarg],"atype") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_atype;
       if (kindflag != NONE && kindflag != ANGLE)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = ANGLE;
 
     } else if (strcmp(arg[iarg],"datom1") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_datom1;
       if (kindflag != NONE && kindflag != DIHEDRAL)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = DIHEDRAL;
     } else if (strcmp(arg[iarg],"datom2") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_datom2;
       if (kindflag != NONE && kindflag != DIHEDRAL)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = DIHEDRAL;
     } else if (strcmp(arg[iarg],"datom3") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_datom3;
       if (kindflag != NONE && kindflag != DIHEDRAL)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = DIHEDRAL;
     } else if (strcmp(arg[iarg],"datom4") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_datom4;
       if (kindflag != NONE && kindflag != DIHEDRAL)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = DIHEDRAL;
     } else if (strcmp(arg[iarg],"dtype") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_dtype;
       if (kindflag != NONE && kindflag != DIHEDRAL)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = DIHEDRAL;
 
     } else if (strcmp(arg[iarg],"iatom1") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_iatom1;
       if (kindflag != NONE && kindflag != IMPROPER)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = IMPROPER;
     } else if (strcmp(arg[iarg],"iatom2") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_iatom2;
       if (kindflag != NONE && kindflag != IMPROPER)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = IMPROPER;
     } else if (strcmp(arg[iarg],"iatom3") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_iatom3;
       if (kindflag != NONE && kindflag != IMPROPER)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = IMPROPER;
     } else if (strcmp(arg[iarg],"iatom4") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_iatom4;
       if (kindflag != NONE && kindflag != IMPROPER)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = IMPROPER;
     } else if (strcmp(arg[iarg],"itype") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_itype;
       if (kindflag != NONE && kindflag != IMPROPER)
         error->all(FLERR,
                    "Compute property/local cannot use these inputs together");
       kindflag = IMPROPER;
 
     } else error->all(FLERR,
                       "Invalid keyword in compute property/local command");
   }
 
   // error check
 
   if (atom->molecular == 2 && (kindflag == BOND || kindflag == ANGLE ||
                                kindflag == DIHEDRAL || kindflag == IMPROPER))
     error->all(FLERR,"Compute property/local does not (yet) work "
                "with atom_style template");
 
   if (kindflag == BOND && atom->avec->bonds_allow == 0)
     error->all(FLERR,
                "Compute property/local for property that isn't allocated");
   if (kindflag == ANGLE && atom->avec->angles_allow == 0)
     error->all(FLERR,
                "Compute property/local for property that isn't allocated");
   if (kindflag == DIHEDRAL && atom->avec->dihedrals_allow == 0)
     error->all(FLERR,
                "Compute property/local for property that isn't allocated");
   if (kindflag == IMPROPER && atom->avec->impropers_allow == 0)
     error->all(FLERR,
                "Compute property/local for property that isn't allocated");
 
   nmax = 0;
   vector = NULL;
   array = NULL;
   indices = NULL;
 }
 
 /* ---------------------------------------------------------------------- */
 
 ComputePropertyLocal::~ComputePropertyLocal()
 {
   delete [] pack_choice;
   memory->destroy(vector);
   memory->destroy(array);
   memory->destroy(indices);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::init()
 {
   if (kindflag == NEIGH || kindflag == PAIR) {
     if (force->pair == NULL)
       error->all(FLERR,"No pair style is defined for compute property/local");
     if (force->pair->single_enable == 0)
       error->all(FLERR,"Pair style does not support compute property/local");
   }
 
   // for NEIGH/PAIR need an occasional half neighbor list
 
   if (kindflag == NEIGH || kindflag == PAIR) {
     int irequest = neighbor->request((void *) this);
     neighbor->requests[irequest]->pair = 0;
     neighbor->requests[irequest]->compute = 1;
     neighbor->requests[irequest]->occasional = 1;
   }
 
   // do initial memory allocation so that memory_usage() is correct
   // cannot be done yet for NEIGH/PAIR, since neigh list does not exist
 
   if (kindflag == NEIGH) ncount = 0;
   else if (kindflag == PAIR) ncount = 0;
   else if (kindflag == BOND) ncount = count_bonds(0);
   else if (kindflag == ANGLE) ncount = count_angles(0);
   else if (kindflag == DIHEDRAL) ncount = count_dihedrals(0);
   else if (kindflag == IMPROPER) ncount = count_impropers(0);
 
   if (ncount > nmax) reallocate(ncount);
   size_local_rows = ncount;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::init_list(int id, NeighList *ptr)
 {
   list = ptr;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::compute_local()
 {
   invoked_local = update->ntimestep;
 
   // count local entries and generate list of indices
 
   if (kindflag == NEIGH) ncount = count_pairs(0,0);
   else if (kindflag == PAIR) ncount = count_pairs(0,1);
   else if (kindflag == BOND) ncount = count_bonds(0);
   else if (kindflag == ANGLE) ncount = count_angles(0);
   else if (kindflag == DIHEDRAL) ncount = count_dihedrals(0);
   else if (kindflag == IMPROPER) ncount = count_impropers(0);
 
   if (ncount > nmax) reallocate(ncount);
   size_local_rows = ncount;
 
   if (kindflag == NEIGH) ncount = count_pairs(1,0);
   else if (kindflag == PAIR) ncount = count_pairs(1,1);
   else if (kindflag == BOND) ncount = count_bonds(1);
   else if (kindflag == ANGLE) ncount = count_angles(1);
   else if (kindflag == DIHEDRAL) ncount = count_dihedrals(1);
   else if (kindflag == IMPROPER) ncount = count_impropers(1);
 
   // fill vector or array with local values
 
   if (nvalues == 1) {
     buf = vector;
     (this->*pack_choice[0])(0);
   } else {
     if (array) buf = &array[0][0];
     for (int n = 0; n < nvalues; n++)
       (this->*pack_choice[n])(n);
   }
 }
 
 /* ----------------------------------------------------------------------
    count pairs and compute pair info on this proc
    only count pair once if newton_pair is off
    both atom I,J must be in group
    if allflag is set, compute requested info about pair
    if forceflag = 1, pair must be within force cutoff, else neighbor cutoff
 ------------------------------------------------------------------------- */
 
 int ComputePropertyLocal::count_pairs(int allflag, int forceflag)
 {
   int i,j,m,ii,jj,inum,jnum,itype,jtype;
+  tagint itag,jtag;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   int *ilist,*jlist,*numneigh,**firstneigh;
 
   double **x = atom->x;
+  tagint *tag = atom->tag;
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
   int newton_pair = force->newton_pair;
 
   // invoke half neighbor list (will copy or build if necessary)
 
   if (allflag == 0) neighbor->build_one(list);
 
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
   // loop over neighbors of my atoms
   // skip if I or J are not in group
+  // for newton = 0 and J = ghost atom,
+  //   need to insure I,J pair is only output by one proc
+  //   use same itag,jtag logic as in Neighbor::neigh_half_nsq()
 
   double **cutsq = force->pair->cutsq;
 
   m = 0;
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
     if (!(mask[i] & groupbit)) continue;
 
     xtmp = x[i][0];
     ytmp = x[i][1];
     ztmp = x[i][2];
+    itag = tag[i];
     itype = type[i];
     jlist = firstneigh[i];
     jnum = numneigh[i];
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
       j &= NEIGHMASK;
 
       if (!(mask[j] & groupbit)) continue;
-      if (newton_pair == 0 && j >= nlocal) continue;
+
+      // itag = jtag is possible for long cutoffs that include images of self
+
+      if (newton_pair == 0 && j >= nlocal) {
+        jtag = tag[j];
+        if (itag > jtag) {
+          if ((itag+jtag) % 2 == 0) continue;
+        } else if (itag < jtag) {
+          if ((itag+jtag) % 2 == 1) continue;
+        } else {
+          if (x[j][2] < ztmp) continue;
+          if (x[j][2] == ztmp) {
+            if (x[j][1] < ytmp) continue;
+            if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
+          }
+        }
+      }
 
       delx = xtmp - x[j][0];
       dely = ytmp - x[j][1];
       delz = ztmp - x[j][2];
       rsq = delx*delx + dely*dely + delz*delz;
       jtype = type[j];
       if (forceflag && rsq >= cutsq[itype][jtype]) continue;
 
       if (allflag) {
         indices[m][0] = i;
         indices[m][1] = j;
       }
       m++;
     }
   }
 
   return m;
 }
 
 /* ----------------------------------------------------------------------
    count bonds on this proc
    only count bond once if newton_bond is off
    all atoms in interaction must be in group
    all atoms in interaction must be known to proc
    if bond is deleted (type = 0), do not count
    if bond is turned off (type < 0), still count
 ------------------------------------------------------------------------- */
 
 int ComputePropertyLocal::count_bonds(int flag)
 {
   int i,atom1,atom2;
 
   int *num_bond = atom->num_bond;
   tagint **bond_atom = atom->bond_atom;
   int **bond_type = atom->bond_type;
   tagint *tag = atom->tag;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
   int newton_bond = force->newton_bond;
 
   int m = 0;
   for (atom1 = 0; atom1 < nlocal; atom1++) {
     if (!(mask[atom1] & groupbit)) continue;
     for (i = 0; i < num_bond[atom1]; i++) {
       atom2 = atom->map(bond_atom[atom1][i]);
       if (atom2 < 0 || !(mask[atom2] & groupbit)) continue;
       if (newton_bond == 0 && tag[atom1] > tag[atom2]) continue;
       if (bond_type[atom1][i] == 0) continue;
 
       if (flag) {
         indices[m][0] = atom1;
         indices[m][1] = i;
       }
       m++;
     }
   }
 
   return m;
 }
 
 /* ----------------------------------------------------------------------
    count angles on this proc
    only count if 2nd atom is the one storing the angle
    all atoms in interaction must be in group
    all atoms in interaction must be known to proc
    if angle is deleted (type = 0), do not count
    if angle is turned off (type < 0), still count
 ------------------------------------------------------------------------- */
 
 int ComputePropertyLocal::count_angles(int flag)
 {
   int i,atom1,atom2,atom3;
 
   int *num_angle = atom->num_angle;
   tagint **angle_atom1 = atom->angle_atom1;
   tagint **angle_atom2 = atom->angle_atom2;
   tagint **angle_atom3 = atom->angle_atom3;
   int **angle_type = atom->angle_type;
   tagint *tag = atom->tag;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
   int m = 0;
   for (atom2 = 0; atom2 < nlocal; atom2++) {
     if (!(mask[atom2] & groupbit)) continue;
     for (i = 0; i < num_angle[atom2]; i++) {
       if (tag[atom2] != angle_atom2[atom2][i]) continue;
       atom1 = atom->map(angle_atom1[atom2][i]);
       if (atom1 < 0 || !(mask[atom1] & groupbit)) continue;
       atom3 = atom->map(angle_atom3[atom2][i]);
       if (atom3 < 0 || !(mask[atom3] & groupbit)) continue;
       if (angle_type[atom2][i] == 0) continue;
 
       if (flag) {
         indices[m][0] = atom2;
         indices[m][1] = i;
       }
       m++;
     }
   }
 
   return m;
 }
 
 /* ----------------------------------------------------------------------
    count dihedrals on this proc
    only count if 2nd atom is the one storing the dihedral
    all atoms in interaction must be in group
    all atoms in interaction must be known to proc
 ------------------------------------------------------------------------- */
 
 int ComputePropertyLocal::count_dihedrals(int flag)
 {
   int i,atom1,atom2,atom3,atom4;
 
   int *num_dihedral = atom->num_dihedral;
   tagint **dihedral_atom1 = atom->dihedral_atom1;
   tagint **dihedral_atom2 = atom->dihedral_atom2;
   tagint **dihedral_atom3 = atom->dihedral_atom3;
   tagint **dihedral_atom4 = atom->dihedral_atom4;
   tagint *tag = atom->tag;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
   int m = 0;
   for (atom2 = 0; atom2 < nlocal; atom2++) {
     if (!(mask[atom2] & groupbit)) continue;
     for (i = 0; i < num_dihedral[atom2]; i++) {
       if (tag[atom2] != dihedral_atom2[atom2][i]) continue;
       atom1 = atom->map(dihedral_atom1[atom2][i]);
       if (atom1 < 0 || !(mask[atom1] & groupbit)) continue;
       atom3 = atom->map(dihedral_atom3[atom2][i]);
       if (atom3 < 0 || !(mask[atom3] & groupbit)) continue;
       atom4 = atom->map(dihedral_atom4[atom2][i]);
       if (atom4 < 0 || !(mask[atom4] & groupbit)) continue;
 
       if (flag) {
         indices[m][0] = atom2;
         indices[m][1] = i;
       }
       m++;
     }
   }
 
   return m;
 }
 
 /* ----------------------------------------------------------------------
    count impropers on this proc
    only count if 2nd atom is the one storing the improper
    all atoms in interaction must be in group
    all atoms in interaction must be known to proc
 ------------------------------------------------------------------------- */
 
 int ComputePropertyLocal::count_impropers(int flag)
 {
   int i,atom1,atom2,atom3,atom4;
 
   int *num_improper = atom->num_improper;
   tagint **improper_atom1 = atom->improper_atom1;
   tagint **improper_atom2 = atom->improper_atom2;
   tagint **improper_atom3 = atom->improper_atom3;
   tagint **improper_atom4 = atom->improper_atom4;
   tagint *tag = atom->tag;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
   int m = 0;
   for (atom2 = 0; atom2 < nlocal; atom2++) {
     if (!(mask[atom2] & groupbit)) continue;
     for (i = 0; i < num_improper[atom2]; i++) {
       if (tag[atom2] != improper_atom2[atom2][i]) continue;
       atom1 = atom->map(improper_atom1[atom2][i]);
       if (atom1 < 0 || !(mask[atom1] & groupbit)) continue;
       atom3 = atom->map(improper_atom3[atom2][i]);
       if (atom3 < 0 || !(mask[atom3] & groupbit)) continue;
       atom4 = atom->map(improper_atom4[atom2][i]);
       if (atom4 < 0 || !(mask[atom4] & groupbit)) continue;
 
       if (flag) {
         indices[m][0] = atom2;
         indices[m][1] = i;
       }
       m++;
     }
   }
 
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::reallocate(int n)
 {
   // grow vector or array and indices array
 
   while (nmax < n) nmax += DELTA;
   if (nvalues == 1) {
     memory->destroy(vector);
     memory->create(vector,nmax,"property/local:vector");
     vector_local = vector;
   } else {
     memory->destroy(array);
     memory->create(array,nmax,nvalues,"property/local:array");
     array_local = array;
   }
 
   memory->destroy(indices);
   memory->create(indices,nmax,2,"property/local:indices");
 }
 
 /* ----------------------------------------------------------------------
    memory usage of local data
 ------------------------------------------------------------------------- */
 
 double ComputePropertyLocal::memory_usage()
 {
   double bytes = nmax*nvalues * sizeof(double);
   bytes += nmax*2 * sizeof(int);
   return bytes;
 }
 
 /* ----------------------------------------------------------------------
    one method for every keyword compute property/local can output
    the atom property is packed into buf starting at n with stride nvalues
    customize a new keyword by adding a method
 ------------------------------------------------------------------------- */
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_patom1(int n)
 {
   int i;
   tagint *tag = atom->tag;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][0];
     buf[n] = tag[i];
     n += nvalues;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_patom2(int n)
 {
   int i;
   tagint *tag = atom->tag;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][1];
     buf[n] = tag[i];
     n += nvalues;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_ptype1(int n)
 {
   int i;
   int *type = atom->type;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][0];
     buf[n] = type[i];
     n += nvalues;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_ptype2(int n)
 {
   int i;
   int *type = atom->type;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][1];
     buf[n] = type[i];
     n += nvalues;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_batom1(int n)
 {
   int i;
   tagint *tag = atom->tag;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][0];
     buf[n] = tag[i];
     n += nvalues;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_batom2(int n)
 {
   int i,j;
   tagint **bond_atom = atom->bond_atom;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][0];
     j = indices[m][1];
     buf[n] = bond_atom[i][j];
     n += nvalues;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_btype(int n)
 {
   int i,j;
   int **bond_type = atom->bond_type;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][0];
     j = indices[m][1];
     buf[n] = bond_type[i][j];
     n += nvalues;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_aatom1(int n)
 {
   int i,j;
   tagint **angle_atom1 = atom->angle_atom1;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][0];
     j = indices[m][1];
     buf[n] = angle_atom1[i][j];
     n += nvalues;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_aatom2(int n)
 {
   int i,j;
   tagint **angle_atom2 = atom->angle_atom2;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][0];
     j = indices[m][1];
     buf[n] = angle_atom2[i][j];
     n += nvalues;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_aatom3(int n)
 {
   int i,j;
   tagint **angle_atom3 = atom->angle_atom3;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][0];
     j = indices[m][1];
     buf[n] = angle_atom3[i][j];
     n += nvalues;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_atype(int n)
 {
   int i,j;
   int **angle_type = atom->angle_type;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][0];
     j = indices[m][1];
     buf[n] = angle_type[i][j];
     n += nvalues;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_datom1(int n)
 {
   int i,j;
   tagint **dihedral_atom1 = atom->dihedral_atom1;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][0];
     j = indices[m][1];
     buf[n] = dihedral_atom1[i][j];
     n += nvalues;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_datom2(int n)
 {
   int i,j;
   tagint **dihedral_atom2 = atom->dihedral_atom2;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][0];
     j = indices[m][1];
     buf[n] = dihedral_atom2[i][j];
     n += nvalues;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_datom3(int n)
 {
   int i,j;
   tagint **dihedral_atom3 = atom->dihedral_atom3;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][0];
     j = indices[m][1];
     buf[n] = dihedral_atom3[i][j];
     n += nvalues;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_datom4(int n)
 {
   int i,j;
   tagint **dihedral_atom4 = atom->dihedral_atom4;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][0];
     j = indices[m][1];
     buf[n] = dihedral_atom4[i][j];
     n += nvalues;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_dtype(int n)
 {
   int i,j;
   int **dihedral_type = atom->dihedral_type;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][0];
     j = indices[m][1];
     buf[n] = dihedral_type[i][j];
     n += nvalues;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_iatom1(int n)
 {
   int i,j;
   tagint **improper_atom1 = atom->improper_atom1;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][0];
     j = indices[m][1];
     buf[n] = improper_atom1[i][j];
     n += nvalues;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_iatom2(int n)
 {
   int i,j;
   tagint **improper_atom2 = atom->improper_atom2;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][0];
     j = indices[m][1];
     buf[n] = improper_atom2[i][j];
     n += nvalues;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_iatom3(int n)
 {
   int i,j;
   tagint **improper_atom3 = atom->improper_atom3;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][0];
     j = indices[m][1];
     buf[n] = improper_atom3[i][j];
     n += nvalues;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_iatom4(int n)
 {
   int i,j;
   tagint **improper_atom4 = atom->improper_atom4;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][0];
     j = indices[m][1];
     buf[n] = improper_atom4[i][j];
     n += nvalues;
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePropertyLocal::pack_itype(int n)
 {
   int i,j;
   int **improper_type = atom->improper_type;
 
   for (int m = 0; m < ncount; m++) {
     i = indices[m][0];
     j = indices[m][1];
     buf[n] = improper_type[i][j];
     n += nvalues;
   }
 }