From f865058700c84dbfc8d1d3a9f30724fcc174748f Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 25 Oct 2012 15:21:33 +0000 Subject: [PATCH 01/12] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8999 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- python/examples/demo.py | 19 +++++++------------ 1 file changed, 7 insertions(+), 12 deletions(-) diff --git a/python/examples/demo.py b/python/examples/demo.py index be38d30167..eb43ec8795 100755 --- a/python/examples/demo.py +++ b/python/examples/demo.py @@ -22,11 +22,6 @@ me = 0 #nprocs = pypar.size() from lammps import lammps -from lammps import LMPINT as INT -from lammps import LMPDOUBLE as DOUBLE -from lammps import LMPIPTR as IPTR -from lammps import LMPDPTR as DPTR -from lammps import LMPDPTRPTR as DPTRPTR lmp = lammps() @@ -36,9 +31,9 @@ lmp.file("in.demo") if me == 0: print "\nPython output:" -natoms = lmp.extract_global("natoms",DOUBLE) -mass = lmp.extract_atom("mass",DPTR) -x = lmp.extract_atom("x",DPTRPTR) +natoms = lmp.extract_global("natoms",0) +mass = lmp.extract_atom("mass",2) +x = lmp.extract_atom("x",3) print "Natoms, mass, x[0][0] coord =",natoms,mass[1],x[0][0] temp = lmp.extract_compute("thermo_temp",0,0) @@ -53,13 +48,13 @@ print "Velocity component from atom-style variable =",vy[1] natoms = lmp.get_natoms() print "Natoms from get_natoms =",natoms -xc = lmp.get_coords() -print "Global coords from get_coords =",xc[0],xc[1],xc[31] +xc = lmp.gather_atoms("x",1,3) +print "Global coords from gather_atoms =",xc[0],xc[1],xc[31] xc[0] = xc[0] + 1.0 -lmp.put_coords(xc) +lmp.scatter_atoms("x",1,3,xc) -print "Changed x[0][0] via put_coords =",x[0][0] +print "Changed x[0][0] via scatter_atoms =",x[0][0] # uncomment if running in parallel via Pypar #print "Proc %d out of %d procs has" % (me,nprocs), lmp From e2283d16092b335d7578951870861cdfc9d32823 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 25 Oct 2012 15:21:37 +0000 Subject: [PATCH 02/12] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9000 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- python/lammps.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/python/lammps.py b/python/lammps.py index 4485bb6b6d..bbd03c593a 100644 --- a/python/lammps.py +++ b/python/lammps.py @@ -88,11 +88,11 @@ class lammps: self.lib.lammps_extract_compute.restype = POINTER(c_double) ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type) return ptr[0] - elif type == 1: + if type == 1: self.lib.lammps_extract_compute.restype = POINTER(c_double) ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type) return ptr - elif type == 2: + if type == 2: self.lib.lammps_extract_compute.restype = POINTER(POINTER(c_double)) ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type) return ptr @@ -109,11 +109,11 @@ class lammps: result = ptr[0] self.lib.lammps_free(ptr) return result - elif type == 1: + if type == 1: self.lib.lammps_extract_fix.restype = POINTER(c_double) ptr = self.lib.lammps_extract_fix(self.lmp,id,style,type,i,j) return ptr - elif type == 2: + if type == 2: self.lib.lammps_extract_fix.restype = POINTER(POINTER(c_double)) ptr = self.lib.lammps_extract_fix(self.lmp,id,style,type,i,j) return ptr From 13dd457ae76f7c158473d137a62e2d20ec3509b6 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 25 Oct 2012 15:23:26 +0000 Subject: [PATCH 03/12] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9001 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- python/lammps.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/lammps.py b/python/lammps.py index bbd03c593a..ac69101168 100644 --- a/python/lammps.py +++ b/python/lammps.py @@ -23,8 +23,8 @@ class lammps: # if name = "g++", load liblammps_g++.so try: - if not name: self.lib = CDLL("liblammps.so") - else: self.lib = CDLL("liblammps_%s.so" % name) + if not name: self.lib = CDLL("liblammps.so",RTLD_GLOBAL) + else: self.lib = CDLL("liblammps_%s.so" % name,RTLD_GLOBAL) except: type,value,tb = sys.exc_info() traceback.print_exception(type,value,tb) From 6dcb5fb23d351bc95a2b5328bbff5e19f4a8b67c Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 25 Oct 2012 15:35:14 +0000 Subject: [PATCH 04/12] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9002 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- examples/COUPLE/README | 2 +- src/Makefile.shlib | 4 +- src/REPLICA/verlet_split.cpp | 18 ++++++- src/USER-OMP/Install.sh | 5 +- src/USER-OMP/Package.sh | 5 +- src/USER-OMP/fix_omp.cpp | 96 +++++++++++++++++++++--------------- src/USER-OMP/thr_omp.cpp | 72 ++------------------------- src/USER-OMP/thr_omp.h | 11 +++-- 8 files changed, 88 insertions(+), 125 deletions(-) diff --git a/examples/COUPLE/README b/examples/COUPLE/README index abd70cb757..137de67601 100644 --- a/examples/COUPLE/README +++ b/examples/COUPLE/README @@ -32,7 +32,7 @@ lammps_spparks grain-growth Monte Carlo with strain via MD, coupling to SPPARKS kinetic MC code library collection of useful inter-code communication routines simple simple example of driver code calling LAMMPS as library -fortran a wrapper on the LAMMPS library API that +fortran a simple wrapper on the LAMMPS library API that can be called from Fortran fortran2 a more sophisticated wrapper on the LAMMPS library API that can be called from Fortran diff --git a/src/Makefile.shlib b/src/Makefile.shlib index 71483a1c87..dd4976d536 100644 --- a/src/Makefile.shlib +++ b/src/Makefile.shlib @@ -7,9 +7,9 @@ SHELL = /bin/sh ROOT = lammps EXE = lib$(ROOT)_$@.so -SRC = angle.cpp angle_charmm.cpp angle_class2.cpp angle_cosine.cpp angle_cosine_delta.cpp angle_cosine_periodic.cpp angle_cosine_squared.cpp angle_harmonic.cpp angle_hybrid.cpp angle_table.cpp atom.cpp atom_map.cpp atom_vec.cpp atom_vec_angle.cpp atom_vec_atomic.cpp atom_vec_bond.cpp atom_vec_charge.cpp atom_vec_dipole.cpp atom_vec_ellipsoid.cpp atom_vec_full.cpp atom_vec_hybrid.cpp atom_vec_line.cpp atom_vec_molecular.cpp atom_vec_peri.cpp atom_vec_sphere.cpp atom_vec_tri.cpp balance.cpp bond.cpp bond_class2.cpp bond_fene.cpp bond_fene_expand.cpp bond_harmonic.cpp bond_hybrid.cpp bond_morse.cpp bond_nonlinear.cpp bond_quartic.cpp bond_table.cpp change_box.cpp comm.cpp compute.cpp compute_angle_local.cpp compute_atom_molecule.cpp compute_bond_local.cpp compute_centro_atom.cpp compute_cluster_atom.cpp compute_cna_atom.cpp compute_com.cpp compute_com_molecule.cpp compute_contact_atom.cpp compute_coord_atom.cpp compute_damage_atom.cpp compute_dihedral_local.cpp compute_displace_atom.cpp compute_erotate_asphere.cpp compute_erotate_sphere.cpp compute_erotate_sphere_atom.cpp compute_event_displace.cpp compute_group_group.cpp compute_gyration.cpp compute_gyration_molecule.cpp compute_heat_flux.cpp compute_improper_local.cpp compute_ke.cpp compute_ke_atom.cpp compute_msd.cpp compute_msd_molecule.cpp compute_pair.cpp compute_pair_local.cpp compute_pe.cpp compute_pe_atom.cpp compute_pressure.cpp compute_property_atom.cpp compute_property_local.cpp compute_property_molecule.cpp compute_rdf.cpp compute_reduce.cpp compute_reduce_region.cpp compute_slice.cpp compute_stress_atom.cpp compute_temp.cpp compute_temp_asphere.cpp compute_temp_com.cpp compute_temp_deform.cpp compute_temp_partial.cpp compute_temp_profile.cpp compute_temp_ramp.cpp compute_temp_region.cpp compute_temp_sphere.cpp compute_ti.cpp create_atoms.cpp create_box.cpp delete_atoms.cpp delete_bonds.cpp dihedral.cpp dihedral_charmm.cpp dihedral_class2.cpp dihedral_harmonic.cpp dihedral_helix.cpp dihedral_hybrid.cpp dihedral_multi_harmonic.cpp dihedral_opls.cpp displace_atoms.cpp domain.cpp dump.cpp dump_atom.cpp dump_cfg.cpp dump_custom.cpp dump_dcd.cpp dump_image.cpp dump_local.cpp dump_xtc.cpp dump_xyz.cpp error.cpp ewald.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_gravity.cpp fix_heat.cpp fix_indent.cpp fix_langevin.cpp fix_lineforce.cpp fix_minimize.cpp fix_momentum.cpp fix_move.cpp fix_msst.cpp fix_neb.cpp fix_nh.cpp fix_nh_asphere.cpp fix_nh_sphere.cpp fix_nph.cpp fix_nph_asphere.cpp fix_nph_sphere.cpp fix_nphug.cpp fix_npt.cpp fix_npt_asphere.cpp fix_npt_sphere.cpp fix_nve.cpp fix_nve_asphere.cpp fix_nve_asphere_noforce.cpp fix_nve_limit.cpp fix_nve_line.cpp fix_nve_noforce.cpp fix_nve_sphere.cpp fix_nve_tri.cpp fix_nvt.cpp fix_nvt_asphere.cpp fix_nvt_sllod.cpp fix_nvt_sphere.cpp fix_orient_fcc.cpp fix_peri_neigh.cpp fix_planeforce.cpp fix_pour.cpp fix_press_berendsen.cpp fix_print.cpp fix_qeq_comb.cpp fix_read_restart.cpp fix_recenter.cpp fix_respa.cpp fix_restrain.cpp fix_rigid.cpp fix_rigid_nve.cpp fix_rigid_nvt.cpp fix_setforce.cpp fix_shake.cpp fix_shear_history.cpp fix_spring.cpp fix_spring_rg.cpp fix_spring_self.cpp fix_srd.cpp fix_store_force.cpp fix_store_state.cpp fix_temp_berendsen.cpp fix_temp_rescale.cpp fix_thermal_conductivity.cpp fix_tmd.cpp fix_ttm.cpp fix_viscosity.cpp fix_viscous.cpp fix_wall.cpp fix_wall_colloid.cpp fix_wall_gran.cpp fix_wall_harmonic.cpp fix_wall_lj126.cpp fix_wall_lj93.cpp fix_wall_piston.cpp fix_wall_reflect.cpp fix_wall_region.cpp fix_wall_srd.cpp force.cpp group.cpp image.cpp improper.cpp improper_class2.cpp improper_cvff.cpp improper_harmonic.cpp improper_hybrid.cpp improper_umbrella.cpp input.cpp integrate.cpp irregular.cpp kspace.cpp lammps.cpp lattice.cpp library.cpp math_extra.cpp memory.cpp min.cpp min_cg.cpp min_fire.cpp min_hftn.cpp min_linesearch.cpp min_quickmin.cpp min_sd.cpp minimize.cpp modify.cpp neb.cpp neigh_bond.cpp neigh_derive.cpp neigh_full.cpp neigh_gran.cpp neigh_half_bin.cpp neigh_half_multi.cpp neigh_half_nsq.cpp neigh_list.cpp neigh_request.cpp neigh_respa.cpp neigh_stencil.cpp neighbor.cpp output.cpp pair.cpp pair_adp.cpp pair_airebo.cpp pair_beck.cpp pair_bop.cpp pair_born.cpp pair_born_coul_long.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_colloid.cpp pair_comb.cpp pair_coul_cut.cpp pair_coul_debye.cpp pair_coul_long.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_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_long.cpp pair_lj_cut_coul_long_opt.cpp pair_lj_cut_coul_long_tip4p.cpp pair_lj_cut_coul_long_tip4p_opt.cpp pair_lj_cut_opt.cpp pair_lj_expand.cpp pair_lj_gromacs.cpp pair_lj_gromacs_coul_gromacs.cpp pair_lj_smooth.cpp pair_lj_smooth_linear.cpp pair_lubricate.cpp pair_lubricateU.cpp pair_lubricateU_poly.cpp pair_lubricate_poly.cpp pair_morse.cpp pair_morse_opt.cpp pair_peri_lps.cpp pair_peri_pmb.cpp pair_rebo.cpp pair_resquared.cpp pair_soft.cpp pair_sw.cpp pair_table.cpp pair_tersoff.cpp pair_tersoff_zbl.cpp pair_tri_lj.cpp pair_yukawa.cpp pair_yukawa_colloid.cpp pppm.cpp pppm_cg.cpp pppm_old.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 run.cpp set.cpp special.cpp tad.cpp temper.cpp thermo.cpp timer.cpp universe.cpp update.cpp variable.cpp velocity.cpp verlet.cpp verlet_split.cpp write_restart.cpp xdr_compat.cpp +SRC = angle.cpp angle_charmm.cpp angle_class2.cpp angle_cosine.cpp angle_cosine_delta.cpp angle_cosine_periodic.cpp angle_cosine_squared.cpp angle_harmonic.cpp angle_hybrid.cpp angle_table.cpp atom.cpp atom_map.cpp atom_vec.cpp atom_vec_angle.cpp atom_vec_atomic.cpp atom_vec_bond.cpp atom_vec_charge.cpp atom_vec_dipole.cpp atom_vec_ellipsoid.cpp atom_vec_full.cpp atom_vec_hybrid.cpp atom_vec_line.cpp atom_vec_molecular.cpp atom_vec_peri.cpp atom_vec_sphere.cpp atom_vec_tri.cpp balance.cpp bond.cpp bond_class2.cpp bond_fene.cpp bond_fene_expand.cpp bond_harmonic.cpp bond_hybrid.cpp bond_morse.cpp bond_nonlinear.cpp bond_quartic.cpp bond_table.cpp change_box.cpp comm.cpp compute.cpp compute_angle_local.cpp compute_atom_molecule.cpp compute_bond_local.cpp compute_centro_atom.cpp compute_cluster_atom.cpp compute_cna_atom.cpp compute_com.cpp compute_com_molecule.cpp compute_contact_atom.cpp compute_coord_atom.cpp compute_damage_atom.cpp compute_dihedral_local.cpp compute_displace_atom.cpp compute_erotate_asphere.cpp compute_erotate_sphere.cpp compute_erotate_sphere_atom.cpp compute_event_displace.cpp compute_group_group.cpp compute_gyration.cpp compute_gyration_molecule.cpp compute_heat_flux.cpp compute_improper_local.cpp compute_ke.cpp compute_ke_atom.cpp compute_msd.cpp compute_msd_molecule.cpp compute_pair.cpp compute_pair_local.cpp compute_pe.cpp compute_pe_atom.cpp compute_pressure.cpp compute_property_atom.cpp compute_property_local.cpp compute_property_molecule.cpp compute_rdf.cpp compute_reduce.cpp compute_reduce_region.cpp compute_slice.cpp compute_stress_atom.cpp compute_temp.cpp compute_temp_asphere.cpp compute_temp_com.cpp compute_temp_deform.cpp compute_temp_partial.cpp compute_temp_profile.cpp compute_temp_ramp.cpp compute_temp_region.cpp compute_temp_sphere.cpp compute_ti.cpp create_atoms.cpp create_box.cpp delete_atoms.cpp delete_bonds.cpp dihedral.cpp dihedral_charmm.cpp dihedral_class2.cpp dihedral_harmonic.cpp dihedral_helix.cpp dihedral_hybrid.cpp dihedral_multi_harmonic.cpp dihedral_opls.cpp displace_atoms.cpp domain.cpp dump.cpp dump_atom.cpp dump_cfg.cpp dump_custom.cpp dump_dcd.cpp dump_image.cpp dump_local.cpp dump_xtc.cpp dump_xyz.cpp error.cpp finish.cpp fix.cpp fix_adapt.cpp fix_addforce.cpp fix_append_atoms.cpp fix_ave_atom.cpp fix_ave_correlate.cpp fix_ave_histo.cpp fix_ave_spatial.cpp fix_ave_time.cpp fix_aveforce.cpp fix_balance.cpp fix_bond_break.cpp fix_bond_create.cpp fix_bond_swap.cpp fix_box_relax.cpp fix_deform.cpp fix_deposit.cpp fix_drag.cpp fix_dt_reset.cpp fix_efield.cpp fix_enforce2d.cpp fix_evaporate.cpp fix_event.cpp fix_event_prd.cpp fix_event_tad.cpp fix_external.cpp fix_freeze.cpp fix_gcmc.cpp fix_gravity.cpp fix_heat.cpp fix_indent.cpp fix_langevin.cpp fix_lineforce.cpp fix_minimize.cpp fix_momentum.cpp fix_move.cpp fix_msst.cpp fix_neb.cpp fix_nh.cpp fix_nh_asphere.cpp fix_nh_sphere.cpp fix_nph.cpp fix_nph_asphere.cpp fix_nph_sphere.cpp fix_nphug.cpp fix_npt.cpp fix_npt_asphere.cpp fix_npt_sphere.cpp fix_nve.cpp fix_nve_asphere.cpp fix_nve_asphere_noforce.cpp fix_nve_limit.cpp fix_nve_line.cpp fix_nve_noforce.cpp fix_nve_sphere.cpp fix_nve_tri.cpp fix_nvt.cpp fix_nvt_asphere.cpp fix_nvt_sllod.cpp fix_nvt_sphere.cpp fix_orient_fcc.cpp fix_peri_neigh.cpp fix_planeforce.cpp fix_pour.cpp fix_press_berendsen.cpp fix_print.cpp fix_qeq_comb.cpp fix_read_restart.cpp fix_recenter.cpp fix_respa.cpp fix_restrain.cpp fix_rigid.cpp fix_rigid_nh.cpp fix_rigid_nph.cpp fix_rigid_npt.cpp fix_rigid_nve.cpp fix_rigid_nvt.cpp fix_setforce.cpp fix_shake.cpp fix_shear_history.cpp fix_spring.cpp fix_spring_rg.cpp fix_spring_self.cpp fix_srd.cpp fix_store_force.cpp fix_store_state.cpp fix_temp_berendsen.cpp fix_temp_rescale.cpp fix_thermal_conductivity.cpp fix_tmd.cpp fix_ttm.cpp fix_viscosity.cpp fix_viscous.cpp fix_wall.cpp fix_wall_colloid.cpp fix_wall_gran.cpp fix_wall_harmonic.cpp fix_wall_lj126.cpp fix_wall_lj93.cpp fix_wall_piston.cpp fix_wall_reflect.cpp fix_wall_region.cpp fix_wall_srd.cpp force.cpp group.cpp image.cpp improper.cpp improper_class2.cpp improper_cvff.cpp improper_harmonic.cpp improper_hybrid.cpp improper_umbrella.cpp input.cpp integrate.cpp irregular.cpp kspace.cpp lammps.cpp lattice.cpp library.cpp math_extra.cpp memory.cpp min.cpp min_cg.cpp min_fire.cpp min_hftn.cpp min_linesearch.cpp min_quickmin.cpp min_sd.cpp minimize.cpp modify.cpp neb.cpp neigh_bond.cpp neigh_derive.cpp neigh_full.cpp neigh_gran.cpp neigh_half_bin.cpp neigh_half_multi.cpp neigh_half_nsq.cpp neigh_list.cpp neigh_request.cpp neigh_respa.cpp neigh_stencil.cpp neighbor.cpp output.cpp pair.cpp pair_adp.cpp pair_airebo.cpp pair_beck.cpp pair_bop.cpp pair_born.cpp pair_born_coul_wolf.cpp pair_brownian.cpp pair_brownian_poly.cpp pair_buck.cpp pair_buck_coul_cut.cpp pair_colloid.cpp pair_comb.cpp pair_coul_cut.cpp pair_coul_debye.cpp pair_coul_dsf.cpp pair_coul_wolf.cpp pair_dipole_cut.cpp pair_dpd.cpp pair_dpd_tstat.cpp pair_dsmc.cpp pair_eam.cpp pair_eam_alloy.cpp pair_eam_alloy_opt.cpp pair_eam_fs.cpp pair_eam_fs_opt.cpp pair_eam_opt.cpp pair_eim.cpp pair_gauss.cpp pair_gayberne.cpp pair_gran_hertz_history.cpp pair_gran_hooke.cpp pair_gran_hooke_history.cpp pair_hbond_dreiding_lj.cpp pair_hbond_dreiding_morse.cpp pair_hybrid.cpp pair_hybrid_overlay.cpp pair_lcbop.cpp pair_line_lj.cpp pair_lj96_cut.cpp pair_lj_charmm_coul_charmm.cpp pair_lj_charmm_coul_charmm_implicit.cpp pair_lj_class2.cpp pair_lj_class2_coul_cut.cpp pair_lj_class2_coul_long.cpp pair_lj_cubic.cpp pair_lj_cut.cpp pair_lj_cut_coul_cut.cpp pair_lj_cut_coul_debye.cpp pair_lj_cut_coul_dsf.cpp pair_lj_cut_opt.cpp pair_lj_expand.cpp pair_lj_gromacs.cpp pair_lj_gromacs_coul_gromacs.cpp pair_lj_smooth.cpp pair_lj_smooth_linear.cpp pair_lubricate.cpp pair_lubricateU.cpp pair_lubricateU_poly.cpp pair_lubricate_poly.cpp pair_morse.cpp pair_morse_opt.cpp pair_peri_lps.cpp pair_peri_pmb.cpp pair_rebo.cpp pair_resquared.cpp pair_soft.cpp pair_sw.cpp pair_table.cpp pair_tersoff.cpp pair_tersoff_zbl.cpp pair_tri_lj.cpp pair_yukawa.cpp pair_yukawa_colloid.cpp prd.cpp procmap.cpp random_mars.cpp random_park.cpp read_data.cpp read_dump.cpp read_restart.cpp reader.cpp reader_native.cpp reader_xyz.cpp region.cpp region_block.cpp region_cone.cpp region_cylinder.cpp region_intersect.cpp region_plane.cpp region_prism.cpp region_sphere.cpp region_union.cpp replicate.cpp rerun.cpp respa.cpp run.cpp set.cpp special.cpp tad.cpp temper.cpp thermo.cpp timer.cpp universe.cpp update.cpp variable.cpp velocity.cpp verlet.cpp verlet_split.cpp write_restart.cpp xdr_compat.cpp -INC = accelerator_cuda.h accelerator_omp.h angle.h angle_charmm.h angle_class2.h angle_cosine.h angle_cosine_delta.h angle_cosine_periodic.h angle_cosine_squared.h angle_harmonic.h angle_hybrid.h angle_table.h atom.h atom_map.h atom_vec.h atom_vec_angle.h atom_vec_atomic.h atom_vec_bond.h atom_vec_charge.h atom_vec_dipole.h atom_vec_ellipsoid.h atom_vec_full.h atom_vec_hybrid.h atom_vec_line.h atom_vec_molecular.h atom_vec_peri.h atom_vec_sphere.h atom_vec_tri.h balance.h bond.h bond_class2.h bond_fene.h bond_fene_expand.h bond_harmonic.h bond_hybrid.h bond_morse.h bond_nonlinear.h bond_quartic.h bond_table.h change_box.h comm.h compute.h compute_angle_local.h compute_atom_molecule.h compute_bond_local.h compute_centro_atom.h compute_cluster_atom.h compute_cna_atom.h compute_com.h compute_com_molecule.h compute_contact_atom.h compute_coord_atom.h compute_damage_atom.h compute_dihedral_local.h compute_displace_atom.h compute_erotate_asphere.h compute_erotate_sphere.h compute_erotate_sphere_atom.h compute_event_displace.h compute_group_group.h compute_gyration.h compute_gyration_molecule.h compute_heat_flux.h compute_improper_local.h compute_ke.h compute_ke_atom.h compute_msd.h compute_msd_molecule.h compute_pair.h compute_pair_local.h compute_pe.h compute_pe_atom.h compute_pressure.h compute_property_atom.h compute_property_local.h compute_property_molecule.h compute_rdf.h compute_reduce.h compute_reduce_region.h compute_slice.h compute_stress_atom.h compute_temp.h compute_temp_asphere.h compute_temp_com.h compute_temp_deform.h compute_temp_partial.h compute_temp_profile.h compute_temp_ramp.h compute_temp_region.h compute_temp_sphere.h compute_ti.h create_atoms.h create_box.h delete_atoms.h delete_bonds.h dihedral.h dihedral_charmm.h dihedral_class2.h dihedral_harmonic.h dihedral_helix.h dihedral_hybrid.h dihedral_multi_harmonic.h dihedral_opls.h displace_atoms.h domain.h dump.h dump_atom.h dump_cfg.h dump_custom.h dump_dcd.h dump_image.h dump_local.h dump_xtc.h dump_xyz.h error.h ewald.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_gravity.h fix_heat.h fix_indent.h fix_langevin.h fix_lineforce.h fix_minimize.h fix_momentum.h fix_move.h fix_msst.h fix_neb.h fix_nh.h fix_nh_asphere.h fix_nh_sphere.h fix_nph.h fix_nph_asphere.h fix_nph_sphere.h fix_nphug.h fix_npt.h fix_npt_asphere.h fix_npt_sphere.h fix_nve.h fix_nve_asphere.h fix_nve_asphere_noforce.h fix_nve_limit.h fix_nve_line.h fix_nve_noforce.h fix_nve_sphere.h fix_nve_tri.h fix_nvt.h fix_nvt_asphere.h fix_nvt_sllod.h fix_nvt_sphere.h fix_orient_fcc.h fix_peri_neigh.h fix_planeforce.h fix_pour.h fix_press_berendsen.h fix_print.h fix_qeq_comb.h fix_read_restart.h fix_recenter.h fix_respa.h fix_restrain.h fix_rigid.h fix_rigid_nve.h fix_rigid_nvt.h fix_setforce.h fix_shake.h fix_shear_history.h fix_spring.h fix_spring_rg.h fix_spring_self.h fix_srd.h fix_store_force.h fix_store_state.h fix_temp_berendsen.h fix_temp_rescale.h fix_thermal_conductivity.h fix_tmd.h fix_ttm.h fix_viscosity.h fix_viscous.h fix_wall.h fix_wall_colloid.h fix_wall_gran.h fix_wall_harmonic.h fix_wall_lj126.h fix_wall_lj93.h fix_wall_piston.h fix_wall_reflect.h fix_wall_region.h fix_wall_srd.h force.h group.h image.h improper.h improper_class2.h improper_cvff.h improper_harmonic.h improper_hybrid.h improper_umbrella.h input.h integrate.h irregular.h kissfft.h kspace.h lammps.h lattice.h library.h lmptype.h lmpwindows.h math_const.h math_extra.h memory.h min.h min_cg.h min_fire.h min_hftn.h min_linesearch.h min_quickmin.h min_sd.h minimize.h modify.h neb.h neigh_bond.h neigh_derive.h neigh_full.h neigh_gran.h neigh_half_bin.h neigh_half_multi.h neigh_half_nsq.h neigh_list.h neigh_request.h neigh_respa.h neighbor.h output.h pack.h pair.h pair_adp.h pair_airebo.h pair_beck.h pair_bop.h pair_born.h pair_born_coul_long.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_colloid.h pair_comb.h pair_coul_cut.h pair_coul_debye.h pair_coul_long.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_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_long.h pair_lj_cut_coul_long_opt.h pair_lj_cut_coul_long_tip4p.h pair_lj_cut_coul_long_tip4p_opt.h pair_lj_cut_opt.h pair_lj_expand.h pair_lj_gromacs.h pair_lj_gromacs_coul_gromacs.h pair_lj_smooth.h pair_lj_smooth_linear.h pair_lubricate.h pair_lubricateU.h pair_lubricateU_poly.h pair_lubricate_poly.h pair_morse.h pair_morse_opt.h pair_peri_lps.h pair_peri_pmb.h pair_rebo.h pair_resquared.h pair_soft.h pair_sw.h pair_table.h pair_tersoff.h pair_tersoff_zbl.h pair_tri_lj.h pair_yukawa.h pair_yukawa_colloid.h pointers.h pppm.h pppm_cg.h pppm_old.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 run.h set.h special.h style_angle.h style_atom.h style_bond.h style_command.h style_compute.h style_dihedral.h style_dump.h style_fix.h style_improper.h style_integrate.h style_kspace.h style_minimize.h style_pair.h style_reader.h style_region.h suffix.h tad.h temper.h thermo.h timer.h universe.h update.h variable.h velocity.h verlet.h verlet_split.h version.h write_restart.h xdr_compat.h +INC = accelerator_cuda.h accelerator_omp.h angle.h angle_charmm.h angle_class2.h angle_cosine.h angle_cosine_delta.h angle_cosine_periodic.h angle_cosine_squared.h angle_harmonic.h angle_hybrid.h angle_table.h atom.h atom_map.h atom_masks.h atom_vec.h atom_vec_angle.h atom_vec_atomic.h atom_vec_bond.h atom_vec_charge.h atom_vec_dipole.h atom_vec_ellipsoid.h atom_vec_full.h atom_vec_hybrid.h atom_vec_line.h atom_vec_molecular.h atom_vec_peri.h atom_vec_sphere.h atom_vec_tri.h balance.h bond.h bond_class2.h bond_fene.h bond_fene_expand.h bond_harmonic.h bond_hybrid.h bond_morse.h bond_nonlinear.h bond_quartic.h bond_table.h change_box.h comm.h compute.h compute_angle_local.h compute_atom_molecule.h compute_bond_local.h compute_centro_atom.h compute_cluster_atom.h compute_cna_atom.h compute_com.h compute_com_molecule.h compute_contact_atom.h compute_coord_atom.h compute_damage_atom.h compute_dihedral_local.h compute_displace_atom.h compute_erotate_asphere.h compute_erotate_sphere.h compute_erotate_sphere_atom.h compute_event_displace.h compute_group_group.h compute_gyration.h compute_gyration_molecule.h compute_heat_flux.h compute_improper_local.h compute_ke.h compute_ke_atom.h compute_msd.h compute_msd_molecule.h compute_pair.h compute_pair_local.h compute_pe.h compute_pe_atom.h compute_pressure.h compute_property_atom.h compute_property_local.h compute_property_molecule.h compute_rdf.h compute_reduce.h compute_reduce_region.h compute_slice.h compute_stress_atom.h compute_temp.h compute_temp_asphere.h compute_temp_com.h compute_temp_deform.h compute_temp_partial.h compute_temp_profile.h compute_temp_ramp.h compute_temp_region.h compute_temp_sphere.h compute_ti.h create_atoms.h create_box.h delete_atoms.h delete_bonds.h dihedral.h dihedral_charmm.h dihedral_class2.h dihedral_harmonic.h dihedral_helix.h dihedral_hybrid.h dihedral_multi_harmonic.h dihedral_opls.h displace_atoms.h domain.h dump.h dump_atom.h dump_cfg.h dump_custom.h dump_dcd.h dump_image.h dump_local.h dump_xtc.h dump_xyz.h error.h finish.h fix.h fix_adapt.h fix_addforce.h fix_append_atoms.h fix_ave_atom.h fix_ave_correlate.h fix_ave_histo.h fix_ave_spatial.h fix_ave_time.h fix_aveforce.h fix_balance.h fix_bond_break.h fix_bond_create.h fix_bond_swap.h fix_box_relax.h fix_deform.h fix_deposit.h fix_drag.h fix_dt_reset.h fix_efield.h fix_enforce2d.h fix_evaporate.h fix_event.h fix_event_prd.h fix_event_tad.h fix_external.h fix_freeze.h fix_gcmc.h fix_gravity.h fix_heat.h fix_indent.h fix_langevin.h fix_lineforce.h fix_minimize.h fix_momentum.h fix_move.h fix_msst.h fix_neb.h fix_nh.h fix_nh_asphere.h fix_nh_sphere.h fix_nph.h fix_nph_asphere.h fix_nph_sphere.h fix_nphug.h fix_npt.h fix_npt_asphere.h fix_npt_sphere.h fix_nve.h fix_nve_asphere.h fix_nve_asphere_noforce.h fix_nve_limit.h fix_nve_line.h fix_nve_noforce.h fix_nve_sphere.h fix_nve_tri.h fix_nvt.h fix_nvt_asphere.h fix_nvt_sllod.h fix_nvt_sphere.h fix_orient_fcc.h fix_peri_neigh.h fix_planeforce.h fix_pour.h fix_press_berendsen.h fix_print.h fix_qeq_comb.h fix_read_restart.h fix_recenter.h fix_respa.h fix_restrain.h fix_rigid.h fix_rigid_nh.h fix_rigid_nph.h fix_rigid_npt.h fix_rigid_nve.h fix_rigid_nvt.h fix_setforce.h fix_shake.h fix_shear_history.h fix_spring.h fix_spring_rg.h fix_spring_self.h fix_srd.h fix_store_force.h fix_store_state.h fix_temp_berendsen.h fix_temp_rescale.h fix_thermal_conductivity.h fix_tmd.h fix_ttm.h fix_viscosity.h fix_viscous.h fix_wall.h fix_wall_colloid.h fix_wall_gran.h fix_wall_harmonic.h fix_wall_lj126.h fix_wall_lj93.h fix_wall_piston.h fix_wall_reflect.h fix_wall_region.h fix_wall_srd.h force.h group.h image.h improper.h improper_class2.h improper_cvff.h improper_harmonic.h improper_hybrid.h improper_umbrella.h input.h integrate.h irregular.h kspace.h lammps.h lattice.h library.h lmptype.h lmpwindows.h math_const.h math_extra.h memory.h min.h min_cg.h min_fire.h min_hftn.h min_linesearch.h min_quickmin.h min_sd.h minimize.h modify.h neb.h neigh_bond.h neigh_derive.h neigh_full.h neigh_gran.h neigh_half_bin.h neigh_half_multi.h neigh_half_nsq.h neigh_list.h neigh_request.h neigh_respa.h neighbor.h output.h pack.h pair.h pair_adp.h pair_airebo.h pair_beck.h pair_bop.h pair_born.h pair_born_coul_wolf.h pair_brownian.h pair_brownian_poly.h pair_buck.h pair_buck_coul_cut.h pair_colloid.h pair_comb.h pair_coul_cut.h pair_coul_debye.h pair_coul_dsf.h pair_coul_wolf.h pair_dipole_cut.h pair_dpd.h pair_dpd_tstat.h pair_dsmc.h pair_eam.h pair_eam_alloy.h pair_eam_alloy_opt.h pair_eam_fs.h pair_eam_fs_opt.h pair_eam_opt.h pair_eim.h pair_gauss.h pair_gayberne.h pair_gran_hertz_history.h pair_gran_hooke.h pair_gran_hooke_history.h pair_hbond_dreiding_lj.h pair_hbond_dreiding_morse.h pair_hybrid.h pair_hybrid_overlay.h pair_lcbop.h pair_line_lj.h pair_lj96_cut.h pair_lj_charmm_coul_charmm.h pair_lj_charmm_coul_charmm_implicit.h pair_lj_class2.h pair_lj_class2_coul_cut.h pair_lj_class2_coul_long.h pair_lj_cubic.h pair_lj_cut.h pair_lj_cut_coul_cut.h pair_lj_cut_coul_debye.h pair_lj_cut_coul_dsf.h pair_lj_cut_opt.h pair_lj_expand.h pair_lj_gromacs.h pair_lj_gromacs_coul_gromacs.h pair_lj_smooth.h pair_lj_smooth_linear.h pair_lubricate.h pair_lubricateU.h pair_lubricateU_poly.h pair_lubricate_poly.h pair_morse.h pair_morse_opt.h pair_peri_lps.h pair_peri_pmb.h pair_rebo.h pair_resquared.h pair_soft.h pair_sw.h pair_table.h pair_tersoff.h pair_tersoff_zbl.h pair_tri_lj.h pair_yukawa.h pair_yukawa_colloid.h pointers.h prd.h procmap.h random_mars.h random_park.h read_data.h read_dump.h read_restart.h reader.h reader_native.h reader_xyz.h region.h region_block.h region_cone.h region_cylinder.h region_intersect.h region_plane.h region_prism.h region_sphere.h region_union.h replicate.h rerun.h respa.h run.h set.h special.h style_angle.h style_atom.h style_bond.h style_command.h style_compute.h style_dihedral.h style_dump.h style_fix.h style_improper.h style_integrate.h style_kspace.h style_minimize.h style_pair.h style_reader.h style_region.h suffix.h tad.h temper.h thermo.h timer.h universe.h update.h variable.h velocity.h verlet.h verlet_split.h version.h write_restart.h xdr_compat.h OBJ = $(SRC:.cpp=.o) diff --git a/src/REPLICA/verlet_split.cpp b/src/REPLICA/verlet_split.cpp index c8754ac681..969768a8bc 100644 --- a/src/REPLICA/verlet_split.cpp +++ b/src/REPLICA/verlet_split.cpp @@ -1,4 +1,4 @@ -/* ------------------------------------------------------------------------- +/* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov @@ -32,6 +32,7 @@ #include "kspace.h" #include "output.h" #include "update.h" +#include "fix.h" #include "modify.h" #include "timer.h" #include "memory.h" @@ -260,7 +261,8 @@ void VerletSplit::setup_minimal(int flag) void VerletSplit::run(int n) { - int nflag,ntimestep,sortflag; + bigint ntimestep; + int nflag,sortflag; // sync both partitions before start timer @@ -272,6 +274,13 @@ void VerletSplit::run(int n) rk_setup(); + // check if OpenMP support fix defined + + Fix *fix_omp; + int ifix = modify->find_fix("package_omp"); + if (ifix < 0) fix_omp = NULL; + else fix_omp = modify->fix[ifix]; + // flags for timestepping iterations int n_post_integrate = modify->n_post_integrate; @@ -360,6 +369,11 @@ void VerletSplit::run(int n) } } else { + + // run FixOMP as sole pre_force fix, if defined + + if (fix_omp) fix_omp->pre_force(vflag); + if (force->kspace) { timer->stamp(); force->kspace->compute(eflag,vflag); diff --git a/src/USER-OMP/Install.sh b/src/USER-OMP/Install.sh index a290fb47bd..8e14d0dc45 100644 --- a/src/USER-OMP/Install.sh +++ b/src/USER-OMP/Install.sh @@ -1,10 +1,9 @@ # Install/unInstall package files in LAMMPS # do not install child files if parent does not exist -for file in *_omp.cpp *_omp.h pppm*proxy.h pppm*proxy.cpp; do +for file in *_omp.cpp *_omp.h ; do # let us see if the "rain man" can count the toothpicks... - ofile=`echo $file | sed -e s,_pppm_tip4p_omp,_long_tip4p_omp, \ - -e s,pppm.\\*_proxy,pppm_omp, -e s,_pppm_omp,_long_omp, \ + ofile=`echo $file | sed \ -e s,\\\\\\(.\\*\\\\\\)_omp\\\\.h,\\\\1.h, \ -e s,\\\\\\(.\\*\\\\\\)_omp\\\\.cpp,\\\\1.cpp,` if (test $1 = 1) then diff --git a/src/USER-OMP/Package.sh b/src/USER-OMP/Package.sh index 6f577b2791..afa02fd4fc 100644 --- a/src/USER-OMP/Package.sh +++ b/src/USER-OMP/Package.sh @@ -4,10 +4,9 @@ # not exist. Do remove OpenMP style files that have no matching # non-OpenMP version installed, e.g. after a package has been # removed -for file in *_omp.cpp *_omp.h pppm*proxy.h pppm*proxy.cpp thr_data.h thr_data.cpp; do +for file in *_omp.cpp *_omp.h thr_data.h thr_data.cpp; do # let us see if the "rain man" can count the toothpicks... - ofile=`echo $file | sed -e s,_pppm_tip4p_omp,_long_tip4p_omp, \ - -e s,pppm.\\*_proxy,pppm_omp, -e s,_pppm_omp,_long_omp, \ + ofile=`echo $file | sed \ -e s,\\\\\\(.\\*\\\\\\)_omp\\\\.\\\\\\(h\\\\\\|cpp\\\\\\),\\\\1.\\\\2,` if (test $file = "thr_omp.h") || (test $file = "thr_omp.cpp") \ || (test $file = "thr_data.h") || (test $file = "thr_data.cpp") then diff --git a/src/USER-OMP/fix_omp.cpp b/src/USER-OMP/fix_omp.cpp index 7e770ce907..ba04a9c2f9 100644 --- a/src/USER-OMP/fix_omp.cpp +++ b/src/USER-OMP/fix_omp.cpp @@ -22,6 +22,7 @@ #include "force.h" #include "neighbor.h" #include "neigh_request.h" +#include "universe.h" #include "update.h" #include "integrate.h" #include "min.h" @@ -194,61 +195,76 @@ void FixOMP::init() if (strstr(update->integrate_style,"respa") != NULL) error->all(FLERR,"Cannot use r-RESPA with /omp styles"); - int check_hybrid; + int check_hybrid, kspace_split; last_pair_hybrid = NULL; last_omp_style = NULL; const char *last_omp_name = NULL; const char *last_hybrid_name = NULL; const char *last_force_name = NULL; + // support for verlet/split operation. + // kspace_split == 0 : regular processing + // kspace_split < 0 : master partition, does not do kspace + // kspace_split > 0 : slave partition, only does kspace + if (strstr(update->integrate_style,"verlet/split") != NULL) { + if (universe->iworld == 0) kspace_split = -1; + else kspace_split = 1; + } else { + kspace_split = 0; + } + // determine which is the last force style with OpenMP // support as this is the one that has to reduce the forces -#define CheckStyleForOMP(name) \ - check_hybrid = 0; \ - if (force->name) { \ - if ( (strcmp(force->name ## _style,"hybrid") == 0) || \ - (strcmp(force->name ## _style,"hybrid/overlay") == 0) ) \ - check_hybrid=1; \ - if (force->name->suffix_flag & Suffix::OMP) { \ - last_force_name = (const char *) #name; \ - last_omp_name = force->name ## _style; \ - last_omp_style = (void *) force->name; \ - } \ +#define CheckStyleForOMP(name) \ + check_hybrid = 0; \ + if (force->name) { \ + if ( (strcmp(force->name ## _style,"hybrid") == 0) || \ + (strcmp(force->name ## _style,"hybrid/overlay") == 0) ) \ + check_hybrid=1; \ + if (force->name->suffix_flag & Suffix::OMP) { \ + last_force_name = (const char *) #name; \ + last_omp_name = force->name ## _style; \ + last_omp_style = (void *) force->name; \ + } \ } -#define CheckHybridForOMP(name,Class) \ - if (check_hybrid) { \ - Class ## Hybrid *style = (Class ## Hybrid *) force->name; \ - for (int i=0; i < style->nstyles; i++) { \ - if (style->styles[i]->suffix_flag & Suffix::OMP) { \ - last_force_name = (const char *) #name; \ - last_omp_name = style->keywords[i]; \ - last_omp_style = style->styles[i]; \ - } \ - } \ +#define CheckHybridForOMP(name,Class) \ + if (check_hybrid) { \ + Class ## Hybrid *style = (Class ## Hybrid *) force->name; \ + for (int i=0; i < style->nstyles; i++) { \ + if (style->styles[i]->suffix_flag & Suffix::OMP) { \ + last_force_name = (const char *) #name; \ + last_omp_name = style->keywords[i]; \ + last_omp_style = style->styles[i]; \ + } \ + } \ } - CheckStyleForOMP(pair); - CheckHybridForOMP(pair,Pair); - if (check_hybrid) { - last_pair_hybrid = last_omp_style; - last_hybrid_name = last_omp_name; + if (kspace_split <= 0) { + CheckStyleForOMP(pair); + CheckHybridForOMP(pair,Pair); + if (check_hybrid) { + last_pair_hybrid = last_omp_style; + last_hybrid_name = last_omp_name; + } + + CheckStyleForOMP(bond); + CheckHybridForOMP(bond,Bond); + + CheckStyleForOMP(angle); + CheckHybridForOMP(angle,Angle); + + CheckStyleForOMP(dihedral); + CheckHybridForOMP(dihedral,Dihedral); + + CheckStyleForOMP(improper); + CheckHybridForOMP(improper,Improper); } - CheckStyleForOMP(bond); - CheckHybridForOMP(bond,Bond); - - CheckStyleForOMP(angle); - CheckHybridForOMP(angle,Angle); - - CheckStyleForOMP(dihedral); - CheckHybridForOMP(dihedral,Dihedral); - - CheckStyleForOMP(improper); - CheckHybridForOMP(improper,Improper); - - CheckStyleForOMP(kspace); + if (kspace_split >= 0) { + CheckStyleForOMP(kspace); + } #undef CheckStyleForOMP #undef CheckHybridForOMP diff --git a/src/USER-OMP/thr_omp.cpp b/src/USER-OMP/thr_omp.cpp index f42b4d327f..ba3088b83b 100644 --- a/src/USER-OMP/thr_omp.cpp +++ b/src/USER-OMP/thr_omp.cpp @@ -147,7 +147,7 @@ void ThrOMP::ev_setup_thr(int eflag, int vflag, int nall, double *eatom, ---------------------------------------------------------------------- */ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, - ThrData *const thr, const int nproxy) + ThrData *const thr) { const int nlocal = lmp->atom->nlocal; const int nghost = lmp->atom->nghost; @@ -172,12 +172,7 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, if (pair->vflag_fdotr) { - if (style == fix->last_pair_hybrid) { - // pair_style hybrid will compute fdotr for us - // but we first need to reduce the forces - data_reduce_thr(&(f[0][0]), nall, nthreads, 3, tid); - need_force_reduce = 0; - } + if (fix->last_pair_hybrid) return; // this is a non-hybrid pair style. compute per thread fdotr if (fix->last_pair_hybrid == NULL) { @@ -205,67 +200,7 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, thr->virial_pair[i] = 0.0; } } - if (eflag & 2) { - data_reduce_thr(&(pair->eatom[0]), nall, nthreads, 1, tid); - } - if (vflag & 4) { - data_reduce_thr(&(pair->vatom[0][0]), nall, nthreads, 6, tid); - } - } - } - break; - case THR_PAIR|THR_PROXY: { - Pair * const pair = lmp->force->pair; - - if (tid >= nproxy && pair->vflag_fdotr) { - - if (fix->last_pair_hybrid) { - if (tid == nproxy) - lmp->error->all(FLERR, - "Cannot use hybrid pair style with kspace proxy"); - else return; - } - - // this is a non-hybrid pair style. compute per thread fdotr - if (fix->last_pair_hybrid == NULL) { - if (lmp->neighbor->includegroup == 0) - thr->virial_fdotr_compute(x, nlocal, nghost, -1); - else - thr->virial_fdotr_compute(x, nlocal, nghost, nfirst); - } - } - - if (evflag) { -#if defined(_OPENMP) -#pragma omp critical -#endif - { - if (tid < nproxy) { - // nothing to collect for kspace proxy threads - // just reset pair accumulators to 0.0. - if (eflag & 1) { - thr->eng_vdwl = 0.0; - thr->eng_coul = 0.0; - } - if (vflag & 3) - for (int i=0; i < 6; ++i) { - thr->virial_pair[i] = 0.0; - } - } else { - if (eflag & 1) { - pair->eng_vdwl += thr->eng_vdwl; - pair->eng_coul += thr->eng_coul; - thr->eng_vdwl = 0.0; - thr->eng_coul = 0.0; - } - if (vflag & 3) - for (int i=0; i < 6; ++i) { - pair->virial[i] += thr->virial_pair[i]; - thr->virial_pair[i] = 0.0; - } - } - } if (eflag & 2) { data_reduce_thr(&(pair->eatom[0]), nall, nthreads, 1, tid); } @@ -439,9 +374,8 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, } break; - case THR_KSPACE|THR_PROXY: // fallthrough case THR_KSPACE: - // nothing to do + // nothing to do. XXX need to add support for per-atom info break; default: diff --git a/src/USER-OMP/thr_omp.h b/src/USER-OMP/thr_omp.h index 19e50031f8..a9f402af97 100644 --- a/src/USER-OMP/thr_omp.h +++ b/src/USER-OMP/thr_omp.h @@ -58,7 +58,8 @@ class ThrOMP { enum {THR_NONE=0,THR_PAIR=1,THR_BOND=1<<1,THR_ANGLE=1<<2, THR_DIHEDRAL=1<<3,THR_IMPROPER=1<<4,THR_KSPACE=1<<5, - THR_CHARMM=1<<6,THR_PROXY=1<<7,THR_HYBRID=1<<8,THR_FIX=1<<9}; + THR_CHARMM=1<<6, /*THR_PROXY=1<<7,*/ THR_HYBRID=1<<8, + THR_FIX=1<<9}; protected: // extra ev_tally setup work for threaded styles @@ -71,7 +72,7 @@ class ThrOMP { // reduce per thread data as needed void reduce_thr(void * const style, const int eflag, const int vflag, - ThrData * const thr, const int nproxy=0); + ThrData * const thr); // thread safe variant error abort support. // signals an error condition in any thread by making @@ -167,14 +168,14 @@ class ThrOMP { // set loop range thread id, and force array offset for threaded runs. static inline void loop_setup_thr(int &ifrom, int &ito, int &tid, - int inum, int nthreads, int nproxy=0) + int inum, int nthreads) { #if defined(_OPENMP) tid = omp_get_thread_num(); // each thread works on a fixed chunk of atoms. - const int idelta = 1 + inum/(nthreads-nproxy); - ifrom = (tid-nproxy)*idelta; + const int idelta = 1 + inum/nthreads; + ifrom = tid*idelta; ito = ((ifrom + idelta) > inum) ? inum : ifrom + idelta; #else tid = 0; From ccbe10ff39d6a1992a0726ce0bcee7817315e862 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 25 Oct 2012 15:42:25 +0000 Subject: [PATCH 05/12] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9003 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- examples/COUPLE/fortran2/LAMMPS.F90 | 2110 ++++++++++++--------------- examples/COUPLE/fortran2/README | 486 +++--- examples/COUPLE/fortran2/in.simple | 9 +- examples/COUPLE/fortran2/simple.f90 | 95 +- 4 files changed, 1292 insertions(+), 1408 deletions(-) diff --git a/examples/COUPLE/fortran2/LAMMPS.F90 b/examples/COUPLE/fortran2/LAMMPS.F90 index 5b8cc6466d..7f19b6b450 100644 --- a/examples/COUPLE/fortran2/LAMMPS.F90 +++ b/examples/COUPLE/fortran2/LAMMPS.F90 @@ -1,1169 +1,941 @@ -!! ----------------------------------------------------------------------- -! LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator -! www.cs.sandia.gov/~sjplimp/lammps.html -! Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories -! -! 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. -!-------------------------------------------------------------------------- - -!! ------------------------------------------------------------------------ -! Contributing author: Karl D. Hammond -! University of Tennessee, Knoxville (USA), 2012 -!-------------------------------------------------------------------------- - -!! LAMMPS, a Fortran 2003 module containing an interface between Fortran -!! programs and the C-style functions in library.cpp that ship with LAMMPS. -!! This file should be accompanied by LAMMPS-wrapper.cpp and LAMMPS-wrapper.h, -!! which define wrapper functions that ease portability and enforce array -!! dimensions. -!! -!! Everything in this module should be 100% portable by way of Fortran 2003's -!! ISO_C_BINDING intrinsic module. See the README for instructions for -!! compilation and use. -!! -!! Here are the PUBLIC functions and subroutines included in this module. -!! subroutine lammps_open (command_line, communicator, ptr) -!! subroutine lammps_open_no_mpi (command_line, ptr) -!! subroutine lammps_close (ptr) -!! subroutine lammps_file (ptr, str) -!! subroutine lammps_command (ptr, str) -!! subroutine lammps_free (ptr) -!! subroutine lammps_extract_global (global, ptr, name) -!! subroutine lammps_extract_atom (atom, ptr, name) -!! subroutine lammps_extract_fix (fix, ptr, id, style, type, i, j) -!! subroutine lammps_extract_compute (compute, ptr, id, style, type) -!! subroutine lammps_extract_variable (variable, ptr, name, group) -!! function lammps_get_natoms (ptr) -!! subroutine lammps_gather_atoms (ptr, name, count, data) -!! subroutine lammps_scatter_atoms (ptr, name, data) - -#define FLERR __FILE__,__LINE__ -! The above line allows for similar error checking as is done with standard -! LAMMPS files. - -module LAMMPS - - use, intrinsic :: ISO_C_binding, only : C_double, C_int, C_ptr, C_char, & - C_NULL_CHAR, C_loc, C_F_pointer, lammps_instance => C_ptr - implicit none - private - public :: lammps_open, lammps_open_no_mpi, lammps_close, lammps_file, & - lammps_command, lammps_free, lammps_extract_global, & - lammps_extract_atom, lammps_extract_compute, lammps_extract_fix, & - lammps_extract_variable, lammps_get_natoms, lammps_gather_atoms, & - lammps_scatter_atoms - public :: lammps_instance - - !! Functions supplemental to the prototypes in library.h. {{{1 - !! The function definitions (in C++) are contained in LAMMPS-wrapper.cpp. - !! I would have written the first in Fortran, but the MPI libraries (which - !! were written in C) have C-based functions to convert from Fortran MPI - !! handles to C MPI handles, and there is no Fortran equivalent for those - !! functions. - interface - subroutine lammps_open_wrapper (argc, argv, communicator, ptr) & - bind (C, name='lammps_open_fortran_wrapper') - import :: C_int, C_ptr - integer (C_int), value :: argc - type (C_ptr), dimension(*) :: argv - integer, value :: communicator - type (C_ptr) :: ptr - end subroutine lammps_open_wrapper - subroutine lammps_actual_error_all (ptr, file, line, str) & - bind (C, name='lammps_error_all') - import :: C_int, C_char, C_ptr - type (C_ptr), value :: ptr - character (kind=C_char), dimension(*), intent(in) :: file, str - integer (C_int), value :: line - end subroutine lammps_actual_error_all - function lammps_get_ntypes (ptr) result (ntypes) & - bind (C, name='lammps_get_ntypes') - import :: C_int, C_ptr - type (C_ptr), value :: ptr - integer (C_int) :: ntypes - end function lammps_get_ntypes - function lammps_actual_extract_compute_vectorsize (ptr, id, style) & - result (vectorsize) bind (C, name='lammps_extract_compute_vectorsize') - import :: C_int, C_char, C_ptr - integer (C_int) :: vectorsize - type (C_ptr), value :: ptr - character (kind=C_char), dimension(*) :: id - integer (C_int), value :: style - end function lammps_actual_extract_compute_vectorsize - subroutine lammps_actual_extract_compute_arraysize (ptr, id, style, & - nrows, ncols) bind (C, name='lammps_extract_compute_arraysize') - import :: C_int, C_char, C_ptr - integer (C_int) :: arraysize - type (C_ptr), value :: ptr - character (kind=C_char), dimension(*) :: id - integer (C_int), value :: style - integer (C_int) :: nrows, ncols - end subroutine lammps_actual_extract_compute_arraysize - function lammps_actual_extract_fix_vectorsize (ptr, id, style) & - result (vectorsize) bind (C, name='lammps_extract_fix_vectorsize') - import :: C_int, C_char, C_ptr - integer (C_int) :: vectorsize - type (C_ptr), value :: ptr - character (kind=C_char), dimension(*) :: id - integer (C_int), value :: style - end function lammps_actual_extract_fix_vectorsize - subroutine lammps_actual_extract_fix_arraysize (ptr, id, style, & - nrows, ncols) bind (C, name='lammps_extract_fix_arraysize') - import :: C_int, C_char, C_ptr - type (C_ptr), value :: ptr - character (kind=C_char), dimension(*) :: id - integer (C_int), value :: style - integer (C_int) :: nrows, ncols - end subroutine lammps_actual_extract_fix_arraysize - end interface - - !! Functions/subroutines defined in library.h and library.cpp {{{1 - interface - subroutine lammps_actual_open_no_mpi (argc, argv, ptr) & - bind (C, name='lammps_open_no_mpi') - import :: C_int, C_ptr - integer (C_int), value :: argc - type (C_ptr), dimension(*) :: argv - type (C_ptr) :: ptr - end subroutine lammps_actual_open_no_mpi - - subroutine lammps_close (ptr) bind (C, name='lammps_close') - import :: C_ptr - type (C_ptr), value :: ptr - end subroutine lammps_close - - subroutine lammps_actual_file (ptr, str) bind (C, name='lammps_file') - import :: C_ptr, C_char - type (C_ptr), value :: ptr - character (kind=C_char), dimension(*) :: str - end subroutine lammps_actual_file - - function lammps_actual_command (ptr, str) result (command) & - bind (C, name='lammps_command') - import :: C_ptr, C_char - type (C_ptr), value :: ptr - character (kind=C_char), dimension(*) :: str - type (C_ptr) :: command - end function lammps_actual_command - - subroutine lammps_free (ptr) bind (C, name='lammps_free') - import :: C_ptr - type (C_ptr), value :: ptr - end subroutine lammps_free - - function lammps_actual_extract_global (ptr, name) & - bind (C, name='lammps_extract_global') result (global) - import :: C_ptr, C_char - type (C_ptr), value :: ptr - character (kind=C_char), dimension(*) :: name - type (C_ptr) :: global - end function lammps_actual_extract_global - - function lammps_actual_extract_atom (ptr, name) & - bind (C, name='lammps_extract_atom') result (atom) - import :: C_ptr, C_char - type (C_ptr), value :: ptr - character (kind=C_char), dimension(*) :: name - type (C_ptr) :: atom - end function lammps_actual_extract_atom - - function lammps_actual_extract_compute (ptr, id, style, type) & - result (compute) bind (C, name='lammps_extract_compute') - import :: C_ptr, C_char, C_int - type (C_ptr), value :: ptr - character (kind=C_char), dimension(*) :: id - integer (C_int), value :: style, type - type (C_ptr) :: compute - end function lammps_actual_extract_compute - - function lammps_actual_extract_fix (ptr, id, style, type, i, j) & - result (fix) bind (C, name='lammps_extract_fix') - import :: C_ptr, C_char, C_int - type (C_ptr), value :: ptr - character (kind=C_char), dimension(*) :: id - integer (C_int), value :: style, type, i, j - type (C_ptr) :: fix - end function lammps_actual_extract_fix - - function lammps_actual_extract_variable (ptr, name, group) & - result (variable) bind (C, name='lammps_extract_variable') - import :: C_ptr, C_char - type (C_ptr), value :: ptr - character (kind=C_char), dimension(*) :: name, group - type (C_ptr) :: variable - end function lammps_actual_extract_variable - - function lammps_get_natoms (ptr) result (natoms) & - bind (C, name='lammps_get_natoms') - import :: C_ptr, C_int - type (C_ptr), value :: ptr - integer (C_int) :: natoms - end function lammps_get_natoms - - subroutine lammps_actual_gather_atoms (ptr, name, type, count, data) & - bind (C, name='lammps_gather_atoms') - import :: C_ptr, C_int, C_char - type (C_ptr), value :: ptr, data - character (kind=C_char), dimension(*) :: name - integer (C_int), value :: type, count - end subroutine lammps_actual_gather_atoms - - subroutine lammps_actual_scatter_atoms (ptr, name, type, count, data) & - bind (C, name='lammps_scatter_atoms') - import :: C_ptr, C_int, C_char - type (C_ptr), value :: ptr, data - character (kind=C_char), dimension(*) :: name - integer (C_int), value :: type, count - end subroutine lammps_actual_scatter_atoms - end interface - - ! Generic functions for the wrappers below {{{1 - - ! Check the dimensions of the arrays these return; they are not always - ! easy to find. Note that I consider returning pointers to arbitrary - ! memory locations with no information as to array size/shape to be - ! extremely sloppy and error-prone. It would appear the Fortran standards - ! committee would agree, as they chose not to allow that sort of nonsense. - - interface lammps_extract_global - module procedure lammps_extract_global_i, lammps_extract_global_r, & - lammps_extract_global_dp - end interface lammps_extract_global - - interface lammps_extract_atom - module procedure lammps_extract_atom_ia, lammps_extract_atom_ra, & - lammps_extract_atom_dpa, lammps_extract_atom_dp2a, & - lammps_extract_atom_r2a - end interface lammps_extract_atom - - interface lammps_extract_compute - module procedure lammps_extract_compute_r, lammps_extract_compute_dp, & - lammps_extract_compute_ra, lammps_extract_compute_dpa, & - lammps_extract_compute_r2a, lammps_extract_compute_dp2a - end interface lammps_extract_compute - - interface lammps_extract_fix - module procedure lammps_extract_fix_r, lammps_extract_fix_dp, & - lammps_extract_fix_ra, lammps_extract_fix_dpa, & - lammps_extract_fix_r2a, lammps_extract_fix_dp2a - end interface lammps_extract_fix - - interface lammps_extract_variable - module procedure lammps_extract_variable_i, & - lammps_extract_variable_dp, & - lammps_extract_variable_r, & - lammps_extract_variable_ra, & - lammps_extract_variable_ia, & - lammps_extract_variable_dpa - end interface lammps_extract_variable - - interface lammps_gather_atoms - module procedure lammps_gather_atoms_ia, lammps_gather_atoms_dpa, & - lammps_gather_atoms_ra - end interface lammps_gather_atoms - - interface lammps_scatter_atoms - module procedure lammps_scatter_atoms_ia, lammps_scatter_atoms_dpa, & - lammps_scatter_atoms_ra - end interface lammps_scatter_atoms - -contains !! Wrapper functions local to this module {{{1 - - subroutine lammps_open (command_line, communicator, ptr) - character (len=*), intent(in) :: command_line - integer, intent(in) :: communicator - type (C_ptr) :: ptr - integer (C_int) :: argc - type (C_ptr), dimension(:), allocatable :: argv - character (kind=C_char), dimension(len_trim(command_line)+1), target :: & - c_command_line - c_command_line = string2Cstring (command_line) - call Cstring2argcargv (c_command_line, argc, argv) - call lammps_open_wrapper (argc, argv, communicator, ptr) - deallocate (argv) - end subroutine lammps_open - -!----------------------------------------------------------------------------- - - subroutine lammps_open_no_mpi (command_line, ptr) - character (len=*), intent(in) :: command_line - type (C_ptr) :: ptr - integer (C_int) :: argc - type (C_ptr), dimension(:), allocatable :: argv - character (kind=C_char), dimension(len_trim(command_line)+1), target :: & - c_command_line - c_command_line = string2Cstring (command_line) - call Cstring2argcargv (c_command_line, argc, argv) - call lammps_actual_open_no_mpi (argc, argv, ptr) - deallocate (argv) - end subroutine lammps_open_no_mpi - -!----------------------------------------------------------------------------- - - subroutine lammps_file (ptr, str) - type (C_ptr) :: ptr - character (len=*) :: str - character (kind=C_char), dimension(len_trim(str)+1) :: Cstr - Cstr = string2Cstring (str) - call lammps_actual_file (ptr, Cstr) - end subroutine lammps_file - -!----------------------------------------------------------------------------- - - subroutine lammps_command (ptr, str) - type (C_ptr) :: ptr - character (len=*) :: str - character (kind=C_char), dimension(len_trim(str)+1) :: Cstr - type (C_ptr) :: dummy - Cstr = string2Cstring (str) - dummy = lammps_actual_command (ptr, Cstr) - end subroutine lammps_command - -!----------------------------------------------------------------------------- - -! lammps_extract_global {{{2 - function lammps_extract_global_Cptr (ptr, name) result (global) - type (C_ptr) :: global - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - character (kind=C_char), dimension(len_trim(name)+1) :: Cname - Cname = string2Cstring (name) - global = lammps_actual_extract_global (ptr, Cname) - end function lammps_extract_global_Cptr - subroutine lammps_extract_global_i (global, ptr, name) - integer, intent(out) :: global - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - type (C_ptr) :: Cptr - integer (C_int), pointer :: Fptr - Cptr = lammps_extract_global_Cptr (ptr, name) - call C_F_pointer (Cptr, Fptr) - global = Fptr - nullify (Fptr) - end subroutine lammps_extract_global_i - subroutine lammps_extract_global_dp (global, ptr, name) - double precision, intent(out) :: global - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - type (C_ptr) :: Cptr - real (C_double), pointer :: Fptr - Cptr = lammps_extract_global_Cptr (ptr, name) - call C_F_pointer (Cptr, Fptr) - global = Fptr - nullify (Fptr) - end subroutine lammps_extract_global_dp - subroutine lammps_extract_global_r (global, ptr, name) - real :: global - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - type (C_ptr) :: Cptr - real (C_double), pointer :: Fptr - Cptr = lammps_extract_global_Cptr (ptr, name) - call C_F_pointer (Cptr, Fptr) - global = real (Fptr) - nullify (Fptr) - end subroutine lammps_extract_global_r - -!----------------------------------------------------------------------------- - -! lammps_extract_atom {{{2 - function lammps_extract_atom_Cptr (ptr, name) result (atom) - type (C_ptr) :: atom - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - character (kind=C_char), dimension(len_trim(name)+1) :: Cname - Cname = string2Cstring (name) - atom = lammps_actual_extract_atom (ptr, Cname) - end function lammps_extract_atom_Cptr - subroutine lammps_extract_atom_ia (atom, ptr, name) - integer, dimension(:), allocatable, intent(out) :: atom - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - type (C_ptr) :: Cptr - integer (C_int), dimension(:), pointer :: Fptr - integer :: nelements - call lammps_extract_global_i (nelements, ptr, 'nlocal') - Cptr = lammps_extract_atom_Cptr (ptr, name) - call C_F_pointer (Cptr, Fptr, (/nelements/)) - if ( .not. associated (Fptr) ) return - allocate (atom(nelements)) - atom = Fptr - nullify (Fptr) - end subroutine lammps_extract_atom_ia - subroutine lammps_extract_atom_dpa (atom, ptr, name) - double precision, dimension(:), allocatable, intent(out) :: atom - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - type (C_ptr) :: Cptr - real (C_double), dimension(:), pointer :: Fptr - integer :: nelements - if ( name == 'mass' ) then - nelements = lammps_get_ntypes (ptr) - else if ( name == 'x' .or. name == 'v' .or. name == 'f' ) then - ! We should not be getting 'x' or 'v' or 'f' here! - call lammps_error_all (ptr, FLERR, 'You cannot extract those atom& - & data (x, v, or f) into a rank 1 array.') - return - else - ! Everything else we can get is probably nlocal units long - call lammps_extract_global_i (nelements, ptr, 'nlocal') - end if - Cptr = lammps_extract_atom_Cptr (ptr, name) - if ( name == 'mass' ) then - call C_F_pointer (Cptr, Fptr, (/nelements + 1/)) - if ( .not. associated (Fptr) ) return - allocate (atom(nelements)) - atom = Fptr(2:) ! LAMMPS starts numbering at 1 (C does not) - else - call C_F_pointer (Cptr, Fptr, (/nelements/)) - if ( .not. associated (Fptr) ) return - allocate (atom(nelements)) - atom = Fptr - end if - nullify (Fptr) - end subroutine lammps_extract_atom_dpa - subroutine lammps_extract_atom_ra (atom, ptr, name) - real, dimension(:), allocatable, intent(out) :: atom - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - double precision, dimension(:), allocatable :: d_atom - call lammps_extract_atom_dpa (d_atom, ptr, name) - allocate (atom(size(d_atom))) - atom = real(d_atom) - deallocate (d_atom) - end subroutine lammps_extract_atom_ra - subroutine lammps_extract_atom_dp2a (atom, ptr, name) - double precision, dimension(:,:), allocatable, intent(out) :: atom - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - type (C_ptr) :: Cptr - integer :: nelements - if ( name /= 'x' .and. name /= 'v' .and. name /= 'f' ) then - call lammps_error_all (ptr, FLERR, 'You cannot extract ' // name // & - ' into a rank 2 array.') - return - end if - Cptr = lammps_extract_atom_Cptr (ptr, name) - call lammps_extract_global_i (nelements, ptr, 'nlocal') - allocate (atom(nelements,3)) - atom = Cdoublestar_to_2darray (Cptr, nelements, 3) - end subroutine lammps_extract_atom_dp2a - subroutine lammps_extract_atom_r2a (atom, ptr, name) - real, dimension(:,:), allocatable, intent(out) :: atom - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - double precision, dimension(:,:), allocatable :: d_atom - call lammps_extract_atom_dp2a (d_atom, ptr, name) - if ( allocated (d_atom) ) then - allocate (atom(size(d_atom,1), size(d_atom,2))) - else - return - end if - atom = real(d_atom) - deallocate (d_atom) - end subroutine lammps_extract_atom_r2a - -!----------------------------------------------------------------------------- - -! lammps_extract_compute {{{2 - function lammps_extract_compute_Cptr (ptr, id, style, type) result (compute) - type (C_ptr) :: compute - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: id - integer, intent(in) :: style, type - integer (kind=C_int) :: Cstyle, Ctype - character (kind=C_char), dimension(len_trim(id)+1) :: Cid - Cid = string2Cstring (id) - Cstyle = style - Ctype = type - compute = lammps_actual_extract_compute (ptr, Cid, Cstyle, Ctype) - end function lammps_extract_compute_Cptr - subroutine lammps_extract_compute_dp (compute, ptr, id, style, type) - double precision, intent(out) :: compute - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: id - integer, intent(in) :: style, type - type (C_ptr) :: Cptr - real (C_double), pointer :: Fptr - ! The only valid values of (style,type) are (0,0) for scalar 'compute' - if ( style /= 0 ) then - call lammps_error_all (ptr, FLERR, 'You cannot pack per-atom/local& - & data into a scalar.') - return - end if - if ( type == 1 ) then - call lammps_error_all (ptr, FLERR, 'You cannot extract a compute& - & vector (rank 1) into a scalar.') - return - else if ( type == 2 ) then - call lammps_error_all (ptr, FLERR, 'You cannot extract a compute& - & array (rank 2) into a scalar.') - return - end if - Cptr = lammps_extract_compute_Cptr (ptr, id, style, type) - call C_F_pointer (Cptr, Fptr) - compute = Fptr - nullify (Fptr) - ! C pointer should not be freed! - end subroutine lammps_extract_compute_dp - subroutine lammps_extract_compute_r (compute, ptr, id, style, type) - real, intent(out) :: compute - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: id - integer, intent(in) :: style, type - double precision :: d_compute - call lammps_extract_compute_dp (d_compute, ptr, id, style, type) - compute = real(d_compute) - end subroutine lammps_extract_compute_r - subroutine lammps_extract_compute_dpa (compute, ptr, id, style, type) - double precision, dimension(:), allocatable, intent(out) :: compute - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: id - integer, intent(in) :: style, type - type (C_ptr) :: Cptr - real (C_double), dimension(:), pointer :: Fptr - integer :: nelements - ! Check for the correct dimensionality - if ( type == 0 ) then - call lammps_error_all (ptr, FLERR, 'You cannot extract a compute& - & scalar (rank 0) into a rank 1 variable.') - return - else if ( type == 2 ) then - call lammps_error_all (ptr, FLERR, 'You cannot extract a compute& - & array (rank 2) into a rank 1 variable.') - return - end if - nelements = lammps_extract_compute_vectorsize (ptr, id, style) - allocate (compute(nelements)) - Cptr = lammps_extract_compute_Cptr (ptr, id, style, type) - call C_F_pointer (Cptr, Fptr, (/nelements/)) - compute = Fptr - nullify (Fptr) - ! C pointer should not be freed - end subroutine lammps_extract_compute_dpa - subroutine lammps_extract_compute_ra (compute, ptr, id, style, type) - real, dimension(:), allocatable, intent(out) :: compute - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: id - integer, intent(in) :: style, type - double precision, dimension(:), allocatable :: d_compute - call lammps_extract_compute_dpa (d_compute, ptr, id, style, type) - allocate (compute(size(d_compute))) - compute = real(d_compute) - deallocate (d_compute) - end subroutine lammps_extract_compute_ra - subroutine lammps_extract_compute_dp2a (compute, ptr, id, style, type) - double precision, dimension(:,:), allocatable, intent(out) :: compute - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: id - integer, intent(in) :: style, type - type (C_ptr) :: Cptr - integer :: nr, nc - ! Check for the correct dimensionality - if ( type == 0 ) then - call lammps_error_all (ptr, FLERR, 'You cannot extract a compute& - & scalar (rank 0) into a rank 2 variable.') - return - else if ( type == 1 ) then - call lammps_error_all (ptr, FLERR, 'You cannot extract a compute& - & array (rank 1) into a rank 2 variable.') - return - end if - call lammps_extract_compute_arraysize (ptr, id, style, nr, nc) - allocate (compute(nr, nc)) - Cptr = lammps_extract_compute_Cptr (ptr, id, style, type) - compute = Cdoublestar_to_2darray (Cptr, nr, nc) - ! C pointer should not be freed - end subroutine lammps_extract_compute_dp2a - subroutine lammps_extract_compute_r2a (compute, ptr, id, style, type) - real, dimension(:,:), allocatable, intent(out) :: compute - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: id - integer, intent(in) :: style, type - double precision, dimension(:,:), allocatable :: d_compute - call lammps_extract_compute_dp2a (d_compute, ptr, id, style, type) - allocate (compute(size(d_compute,1), size(d_compute,2))) - compute = real(d_compute) - deallocate (d_compute) - end subroutine lammps_extract_compute_r2a - -!----------------------------------------------------------------------------- - -! lammps_extract_fix {{{2 - function lammps_extract_fix_Cptr (ptr, id, style, type, i, j) & - result (fix) - type (C_ptr) :: fix - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: id - integer, intent(in) :: style, type, i, j - character (kind=C_char), dimension(len_trim(id)+1) :: Cid - integer (kind=C_int) :: Cstyle, Ctype, Ci, Cj - Cid = string2Cstring (id) - Cstyle = style - Ctype = type - Ci = i - 1 ! This is for consistency with the values from f_ID[i], - Cj = j - 1 ! which is different from what library.cpp uses! - if ( (type >= 1 .and. Ci < 0) .or. & - (type == 2 .and. (Ci < 0 .or. Cj < 0) ) ) then - call lammps_error_all (ptr, FLERR, 'Index out of range in& - & lammps_extract_fix') - end if - fix = lammps_actual_extract_fix (ptr, Cid, Cstyle, Ctype, Ci, Cj) - end function lammps_extract_fix_Cptr - subroutine lammps_extract_fix_dp (fix, ptr, id, style, type, i, j) - double precision, intent(out) :: fix - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: id - integer, intent(in) :: style, type, i, j - type (C_ptr) :: Cptr - real (C_double), pointer :: Fptr - ! Check for the correct dimensionality - if ( style /= 0 ) then - select case (type) - case (0) - call lammps_error_all (ptr, FLERR, 'There is no per-atom or local& - & scalar data available from fixes.') - case (1) - call lammps_error_all (ptr, FLERR, 'You cannot extract a fix''s & - &per-atom/local vector (rank 1) into a scalar.') - case (2) - call lammps_error_all (ptr, FLERR, 'You cannot extract a fix''s & - &per-atom/local array (rank 2) into a scalar.') - case default - call lammps_error_all (ptr, FLERR, 'Invalid extract_fix style& - & value.') - end select - return - end if - Cptr = lammps_extract_fix_Cptr (ptr, id, style, type, i, j) - call C_F_pointer (Cptr, Fptr) - fix = Fptr - nullify (Fptr) - ! Memory is only allocated for "global" fix variables - if ( style == 0 ) call lammps_free (Cptr) - end subroutine lammps_extract_fix_dp - subroutine lammps_extract_fix_r (fix, ptr, id, style, type, i, j) - real, intent(out) :: fix - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: id - integer, intent(in) :: style, type, i, j - double precision :: d_fix - call lammps_extract_fix_dp (d_fix, ptr, id, style, type, i, j) - fix = real(d_fix) - end subroutine lammps_extract_fix_r - subroutine lammps_extract_fix_dpa (fix, ptr, id, style, type, i, j) - double precision, dimension(:), allocatable, intent(out) :: fix - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: id - integer, intent(in) :: style, type, i, j - type (C_ptr) :: Cptr - real (C_double), dimension(:), pointer :: Fptr - integer :: fix_len - ! Check for the correct dimensionality - if ( style == 0 ) then - call lammps_error_all (ptr, FLERR, 'You can''t extract the& - & whole vector from global fix data') - return - else if ( type == 0 ) then - call lammps_error_all (ptr, FLERR, 'You can''t extract a fix& - & scalar into a rank 1 variable') - return - else if ( type == 2 ) then - call lammps_error_all (ptr, FLERR, 'You cannot extract a fix& - & array into a rank 1 variable.') - return - else if ( type /= 1 ) then - call lammps_error_all (ptr, FLERR, 'Invalid type for fix extraction.') - return - end if - fix_len = lammps_extract_fix_vectorsize (ptr, id, style) - allocate (fix(fix_len)) - Cptr = lammps_extract_fix_Cptr (ptr, id, style, type, i, j) - call C_F_pointer (Cptr, Fptr, (/fix_len/)) - fix = Fptr - nullify (Fptr) - ! Memory is only allocated for "global" fix variables - if ( style == 0 ) call lammps_free (Cptr) - end subroutine lammps_extract_fix_dpa - subroutine lammps_extract_fix_ra (fix, ptr, id, style, type, i, j) - real, dimension(:), allocatable, intent(out) :: fix - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: id - integer, intent(in) :: style, type, i, j - double precision, dimension(:), allocatable :: d_fix - call lammps_extract_fix_dpa (d_fix, ptr, id, style, type, i, j) - allocate (fix(size(d_fix))) - fix = real(d_fix) - deallocate (d_fix) - end subroutine lammps_extract_fix_ra - subroutine lammps_extract_fix_dp2a (fix, ptr, id, style, type, i, j) - double precision, dimension(:,:), allocatable, intent(out) :: fix - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: id - integer, intent(in) :: style, type, i, j - type (C_ptr) :: Cptr - integer :: nr, nc - ! Check for the correct dimensionality - if ( style == 0 ) then - call lammps_error_all (ptr, FLERR, 'It is not possible to extract the& - & entire array from global fix data.') - return - else if ( type == 0 ) then - call lammps_error_all (ptr, FLERR, 'You cannot extract a fix& - & scalar (rank 0) into a rank 2 variable.') - return - else if ( type == 1 ) then - call lammps_error_all (ptr, FLERR, 'You cannot extract a fix& - & vector (rank 1) into a rank 2 variable.') - return - end if - call lammps_extract_fix_arraysize (ptr, id, style, nr, nc) - allocate (fix(nr, nc)) - Cptr = lammps_extract_fix_Cptr (ptr, id, style, type, i, j) - fix = Cdoublestar_to_2darray (Cptr, nr, nc) - ! C pointer should not be freed - end subroutine lammps_extract_fix_dp2a - subroutine lammps_extract_fix_r2a (fix, ptr, id, style, type, i, j) - real, dimension(:,:), allocatable, intent(out) :: fix - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: id - integer, intent(in) :: style, type, i, j - double precision, dimension(:,:), allocatable :: d_fix - call lammps_extract_fix_dp2a (d_fix, ptr, id, style, type, i, j) - allocate (fix(size(d_fix,1), size(d_fix,2))) - fix = real(d_fix) - deallocate (d_fix) - end subroutine lammps_extract_fix_r2a - -!----------------------------------------------------------------------------- - -! lammps_extract_variable {{{2 - function lammps_extract_variable_Cptr (ptr, name, group) result (variable) - type (C_ptr) :: ptr, variable - character (len=*) :: name - character (len=*), optional :: group - character (kind=C_char), dimension(len_trim(name)+1) :: Cname - character (kind=C_char), dimension(:), allocatable :: Cgroup - Cname = string2Cstring (name) - if ( present(group) ) then - allocate (Cgroup(len_trim(group)+1)) - Cgroup = string2Cstring (group) - else - allocate (Cgroup(1)) - Cgroup(1) = C_NULL_CHAR - end if - variable = lammps_actual_extract_variable (ptr, Cname, Cgroup) - deallocate (Cgroup) - end function lammps_extract_variable_Cptr - subroutine lammps_extract_variable_i (variable, ptr, name, group) - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - character (len=*), intent(in), optional :: group - integer, intent(out) :: variable - double precision :: d_var - if ( present (group) ) then - call lammps_extract_variable_dp (d_var, ptr, name, group) - else - call lammps_extract_variable_dp (d_var, ptr, name) - end if - variable = nint(d_var) - end subroutine lammps_extract_variable_i - subroutine lammps_extract_variable_dp (variable, ptr, name, group) - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - character (len=*), intent(in), optional :: group - double precision, intent(out) :: variable - type (C_ptr) :: Cptr - real (C_double), pointer :: Fptr - if ( present(group) ) then - Cptr = lammps_extract_variable_Cptr (ptr, name, group) - else - Cptr = lammps_extract_variable_Cptr (ptr, name) - end if - call C_F_pointer (Cptr, Fptr) - variable = Fptr - nullify (Fptr) - call lammps_free (Cptr) - end subroutine lammps_extract_variable_dp - subroutine lammps_extract_variable_r (variable, ptr, name, group) - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - character (len=*), intent(in), optional :: group - real, intent(out) :: variable - double precision :: d_var - if ( present (group) ) then - call lammps_extract_variable_dp (d_var, ptr, name, group) - else - call lammps_extract_variable_dp (d_var, ptr, name) - end if - variable = real(d_var) - end subroutine lammps_extract_variable_r - - subroutine lammps_extract_variable_ia (variable, ptr, name, group) - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - character (len=*), intent(in), optional :: group - integer, dimension(:), allocatable, intent(out) :: variable - double precision, dimension(:), allocatable :: d_var - if ( present (group) ) then - call lammps_extract_variable_dpa (d_var, ptr, name, group) - else - call lammps_extract_variable_dpa (d_var, ptr, name) - end if - allocate (variable(size(d_var))) - variable = nint(d_var) - deallocate (d_var) - end subroutine lammps_extract_variable_ia - subroutine lammps_extract_variable_dpa (variable, ptr, name, group) - double precision, dimension(:), allocatable, intent(out) :: variable - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - character (len=*), intent(in), optional :: group - type (C_ptr) :: Cptr - real (C_double), dimension(:), pointer :: Fptr - integer :: natoms - if ( present(group) ) then - Cptr = lammps_extract_variable_Cptr (ptr, name, group) - else - Cptr = lammps_extract_variable_Cptr (ptr, name) - end if - natoms = lammps_get_natoms (ptr) - allocate (variable(natoms)) - call C_F_pointer (Cptr, Fptr, (/natoms/)) - variable = Fptr - nullify (Fptr) - call lammps_free (Cptr) - end subroutine lammps_extract_variable_dpa - subroutine lammps_extract_variable_ra (variable, ptr, name, group) - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - character (len=*), intent(in), optional :: group - real, dimension(:), allocatable, intent(out) :: variable - double precision, dimension(:), allocatable :: d_var - if ( present (group) ) then - call lammps_extract_variable_dpa (d_var, ptr, name, group) - else - call lammps_extract_variable_dpa (d_var, ptr, name) - end if - allocate (variable(size(d_var))) - variable = real(d_var) - deallocate (d_var) - end subroutine lammps_extract_variable_ra - -!-------------------------------------------------------------------------2}}} - - subroutine lammps_gather_atoms_ia (ptr, name, count, data) - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - integer, intent(in) :: count - integer, dimension(:), allocatable, intent(out) :: data - type (C_ptr) :: Cdata - integer (C_int), dimension(:), pointer :: Fdata - integer (C_int) :: natoms - character (kind=C_char), dimension(len_trim(name)+1) :: Cname - integer (C_int), parameter :: Ctype = 0_C_int - integer (C_int) :: Ccount - natoms = lammps_get_natoms (ptr) - Cname = string2Cstring (name) - if ( count /= 1 .and. count /= 3 ) then - call lammps_error_all (ptr, FLERR, 'lammps_gather_atoms requires& - & count to be either 1 or 3') - else - Ccount = count - end if - allocate ( Fdata(count*natoms) ) - allocate ( data(count*natoms) ) - Cdata = C_loc (Fdata(1)) - call lammps_actual_gather_atoms (ptr, Cname, Ctype, Ccount, Cdata) - data = Fdata - deallocate (Fdata) - end subroutine lammps_gather_atoms_ia - subroutine lammps_gather_atoms_dpa (ptr, name, count, data) - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - integer, intent(in) :: count - double precision, dimension(:), allocatable, intent(out) :: data - type (C_ptr) :: Cdata - real (C_double), dimension(:), pointer :: Fdata - integer (C_int) :: natoms - character (kind=C_char), dimension(len_trim(name)+1) :: Cname - integer (C_int), parameter :: Ctype = 1_C_int - integer (C_int) :: Ccount - natoms = lammps_get_natoms (ptr) - Cname = string2Cstring (name) - if ( count /= 1 .and. count /= 3 ) then - call lammps_error_all (ptr, FLERR, 'lammps_gather_atoms requires& - & count to be either 1 or 3') - else - Ccount = count - end if - allocate ( Fdata(count*natoms) ) - allocate ( data(count*natoms) ) - Cdata = C_loc (Fdata(1)) - call lammps_actual_gather_atoms (ptr, Cname, Ctype, Ccount, Cdata) - data = Fdata(:) - deallocate (Fdata) - end subroutine lammps_gather_atoms_dpa - subroutine lammps_gather_atoms_ra (ptr, name, count, data) - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - integer, intent(in) :: count - real, dimension(:), allocatable, intent(out) :: data - double precision, dimension(:), allocatable :: d_data - call lammps_gather_atoms_dpa (ptr, name, count, d_data) - allocate (data(size(d_data))) - data = d_data - deallocate (d_data) - end subroutine lammps_gather_atoms_ra - -!----------------------------------------------------------------------------- - - subroutine lammps_scatter_atoms_ia (ptr, name, data) - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - integer, dimension(:), intent(in) :: data - integer (kind=C_int) :: natoms, Ccount - integer (kind=C_int), parameter :: Ctype = 0_C_int - character (kind=C_char), dimension(len_trim(name)+1) :: Cname - integer (C_int), dimension(size(data)), target :: Fdata - type (C_ptr) :: Cdata - natoms = lammps_get_natoms (ptr) - Cname = string2Cstring (name) - Ccount = size(data) / natoms - if ( Ccount /= 1 .and. Ccount /= 3 ) & - call lammps_error_all (ptr, FLERR, 'lammps_gather_atoms requires& - & count to be either 1 or 3') - Fdata = data - Cdata = C_loc (Fdata(1)) - call lammps_actual_scatter_atoms (ptr, Cname, Ctype, Ccount, Cdata) - end subroutine lammps_scatter_atoms_ia - subroutine lammps_scatter_atoms_dpa (ptr, name, data) - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - double precision, dimension(:), intent(in) :: data - integer (kind=C_int) :: natoms, Ccount - integer (kind=C_int), parameter :: Ctype = 1_C_int - character (kind=C_char), dimension(len_trim(name)+1) :: Cname - real (C_double), dimension(size(data)), target :: Fdata - type (C_ptr) :: Cdata - natoms = lammps_get_natoms (ptr) - Cname = string2Cstring (name) - Ccount = size(data) / natoms - if ( Ccount /= 1 .and. Ccount /= 3 ) & - call lammps_error_all (ptr, FLERR, 'lammps_gather_atoms requires& - & count to be either 1 or 3') - Fdata = data - Cdata = C_loc (Fdata(1)) - call lammps_actual_scatter_atoms (ptr, Cname, Ctype, Ccount, Cdata) - end subroutine lammps_scatter_atoms_dpa - subroutine lammps_scatter_atoms_ra (ptr, name, data) - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: name - real, dimension(:), intent(in) :: data - double precision, dimension(size(data)) :: d_data - d_data = real (data, kind(d_data)) - call lammps_scatter_atoms_dpa (ptr, name, d_data) - end subroutine lammps_scatter_atoms_ra - -!----------------------------------------------------------------------------- - - function lammps_extract_compute_vectorsize (ptr, id, style) & - result (vectorsize) - integer :: vectorsize - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: id - integer, intent(in) :: style - integer (C_int) :: Cvectorsize, Cstyle - character (kind=C_char), dimension(len_trim(id)+1) :: Cid - Cid = string2Cstring (id) - Cstyle = int(style, C_int) - Cvectorsize = lammps_actual_extract_compute_vectorsize (ptr, Cid, Cstyle) - vectorsize = int(Cvectorsize, kind(vectorsize)) - end function lammps_extract_compute_vectorsize - -!----------------------------------------------------------------------------- - - function lammps_extract_fix_vectorsize (ptr, id, style) & - result (vectorsize) - integer :: vectorsize - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: id - integer, intent(in) :: style - integer (C_int) :: Cvectorsize, Cstyle - character (kind=C_char), dimension(len_trim(id)+1) :: Cid - Cid = string2Cstring (id) - Cstyle = int(style, C_int) - Cvectorsize = lammps_actual_extract_fix_vectorsize (ptr, Cid, Cstyle) - vectorsize = int(Cvectorsize, kind(vectorsize)) - end function lammps_extract_fix_vectorsize - -!----------------------------------------------------------------------------- - - subroutine lammps_extract_compute_arraysize (ptr, id, style, nrows, ncols) - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: id - integer, intent(in) :: style - integer, intent(out) :: nrows, ncols - integer (C_int) :: Cstyle, Cnrows, Cncols - character (kind=C_char), dimension(len_trim(id)+1) :: Cid - Cid = string2Cstring (id) - Cstyle = int (style, C_int) - call lammps_actual_extract_compute_arraysize (ptr, Cid, Cstyle, & - Cnrows, Cncols) - nrows = int (Cnrows, kind(nrows)) - ncols = int (Cncols, kind(ncols)) - end subroutine lammps_extract_compute_arraysize - -!----------------------------------------------------------------------------- - - subroutine lammps_extract_fix_arraysize (ptr, id, style, nrows, ncols) - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: id - integer, intent(in) :: style - integer, intent(out) :: nrows, ncols - integer (C_int) :: Cstyle, Cnrows, Cncols - character (kind=C_char), dimension(len_trim(id)+1) :: Cid - Cid = string2Cstring (id) - Cstyle = int (style, kind(Cstyle)) - call lammps_actual_extract_fix_arraysize (ptr, Cid, Cstyle, & - Cnrows, Cncols) - nrows = int (Cnrows, kind(nrows)) - ncols = int (Cncols, kind(ncols)) - end subroutine lammps_extract_fix_arraysize - -!----------------------------------------------------------------------------- - - subroutine lammps_error_all (ptr, file, line, str) - type (C_ptr), intent(in) :: ptr - character (len=*), intent(in) :: file, str - integer, intent(in) :: line - character (kind=C_char), dimension(len_trim(file)+1) :: Cfile - character (kind=C_char), dimension(len_trim(str)+1) :: Cstr - integer (C_int) :: Cline - Cline = int(line, kind(Cline)) - Cfile = string2Cstring (file) - Cstr = string2Cstring (str) - call lammps_actual_error_all (ptr, Cfile, Cline, Cstr) - end subroutine lammps_error_all - -!----------------------------------------------------------------------------- - -! Locally defined helper functions {{{1 - - pure function string2Cstring (string) result (C_string) - use, intrinsic :: ISO_C_binding, only : C_char, C_NULL_CHAR - character (len=*), intent(in) :: string - character (len=1, kind=C_char) :: C_string (len_trim(string)+1) - integer :: i, n - n = len_trim (string) - forall (i = 1:n) - C_string(i) = string(i:i) - end forall - C_string(n+1) = C_NULL_CHAR - end function string2Cstring - -!----------------------------------------------------------------------------- - - subroutine Cstring2argcargv (Cstring, argc, argv) - !! Converts a C-style string to argc and argv, that is, words in Cstring - !! become C-style strings in argv. IMPORTANT: Cstring is modified by - !! this routine! I would make Cstring local TO this routine and accept - !! a Fortran-style string instead, but we run into scoping and - !! allocation problems that way. This routine assumes the string is - !! null-terminated, as all C-style strings must be. - - character (kind=C_char), dimension(*), target, intent(inout) :: Cstring - integer (C_int), intent(out) :: argc - type (C_ptr), dimension(:), allocatable, intent(out) :: argv - - integer :: StringStart, SpaceIndex, strlen, argnum - - argc = 1_C_int - - ! Find the length of the string - strlen = 1 - do while ( Cstring(strlen) /= C_NULL_CHAR ) - strlen = strlen + 1 - end do - - ! Find the number of non-escaped spaces - SpaceIndex = 2 - do while ( SpaceIndex < strlen ) - if ( Cstring(SpaceIndex) == ' ' .and. & - Cstring(SpaceIndex-1) /= '\' ) then - argc = argc + 1_C_int - ! Find the next non-space character - do while ( Cstring(SpaceIndex+1) == ' ') - SpaceIndex = SpaceIndex + 1 - end do - end if - SpaceIndex = SpaceIndex + 1 - end do - - ! Now allocate memory for argv - allocate (argv(argc)) - - ! Now find the string starting and ending locations - StringStart = 1 - SpaceIndex = 2 - argnum = 1 - do while ( SpaceIndex < strlen ) - if ( Cstring(SpaceIndex) == ' ' .and. & - Cstring(SpaceIndex-1) /= '\' ) then - ! Found a real space => split strings and store this one - Cstring(Spaceindex) = C_NULL_CHAR ! Replaces space with NULL - argv(argnum) = C_loc(Cstring(StringStart)) - argnum = argnum + 1 - ! Find the next non-space character - do while ( Cstring(SpaceIndex+1) == ' ') - SpaceIndex = SpaceIndex + 1 - end do - StringStart = SpaceIndex + 1 - else if ( Cstring(SpaceIndex) == ' ' .and. & - Cstring(SpaceIndex-1) == '\' ) then - ! Escaped space => remove backslash and move rest of array - Cstring(SpaceIndex-1:strlen-1) = Cstring(SpaceIndex:strlen) - strlen = strlen - 1 ! Last character is still C_NULL_CHAR - end if - SpaceIndex = SpaceIndex + 1 - end do - ! Now handle the last argument - argv(argnum) = C_loc(Cstring(StringStart)) - - end subroutine Cstring2argcargv - -!----------------------------------------------------------------------------- - - function Cdoublestar_to_2darray (Carray, nrows, ncolumns) result (Farray) - - ! Take a C/C++ array of pointers to pointers to doubles (sort of like a - ! two-dimensional array, and handled the same way from the programmer's - ! perspective) into a Fortran-style array. Note that columns in C still - ! correspond to columns in Fortran here and the same for rows. - - type (C_ptr), intent(in) :: Carray - integer, intent(in) :: nrows, ncolumns - double precision, dimension(nrows, ncolumns) :: Farray - type (C_ptr), dimension(:), pointer :: C_rows - real (C_double), dimension(:), pointer :: F_row - integer :: i - - ! Convert each "C row pointer" into an array of rows - call C_F_pointer (Carray, C_rows, (/nrows/)) - do i = 1, nrows - ! Convert each C pointer (an entire row) into a Fortran pointer - call C_F_pointer (C_rows(i), F_row, (/ncolumns/)) - Farray (i,:) = real(F_row, kind(0.0D0)) - end do - - end function Cdoublestar_to_2darray -! 1}}} - -end module LAMMPS - -! vim: foldmethod=marker tabstop=3 softtabstop=3 shiftwidth=3 expandtab +!! ----------------------------------------------------------------------- +! LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator +! www.cs.sandia.gov/~sjplimp/lammps.html +! Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories +! +! 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. +!-------------------------------------------------------------------------- + +!! ------------------------------------------------------------------------ +! Contributing author: Karl D. Hammond +! University of Tennessee, Knoxville (USA), 2012 +!-------------------------------------------------------------------------- + +!! LAMMPS, a Fortran 2003 module containing an interface between Fortran +!! programs and the C-style functions in library.cpp that ship with LAMMPS. +!! This file should be accompanied by LAMMPS-wrapper.cpp and LAMMPS-wrapper.h, +!! which define wrapper functions that ease portability and enforce array +!! dimensions. +!! +!! Everything in this module should be 100% portable by way of Fortran 2003's +!! ISO_C_BINDING intrinsic module. See the README for instructions for +!! compilation and use. +!! +!! Here are the PUBLIC functions and subroutines included in this module. +!! subroutine lammps_open (command_line, communicator, ptr) +!! subroutine lammps_open_no_mpi (command_line, ptr) +!! subroutine lammps_close (ptr) +!! subroutine lammps_file (ptr, str) +!! subroutine lammps_command (ptr, str) +!! subroutine lammps_free (ptr) +!! subroutine lammps_extract_global (global, ptr, name) +!! subroutine lammps_extract_atom (atom, ptr, name) +!! subroutine lammps_extract_fix (fix, ptr, id, style, type, i, j) +!! subroutine lammps_extract_compute (compute, ptr, id, style, type) +!! subroutine lammps_extract_variable (variable, ptr, name, group) +!! function lammps_get_natoms (ptr) +!! subroutine lammps_gather_atoms (ptr, name, count, data) +!! subroutine lammps_scatter_atoms (ptr, name, data) + +#define FLERR __FILE__,__LINE__ +! The above line allows for similar error checking as is done with standard +! LAMMPS files. + +module LAMMPS + + use, intrinsic :: ISO_C_binding, only : C_double, C_int, C_ptr, C_char, & + C_NULL_CHAR, C_loc, C_F_pointer, lammps_instance => C_ptr + implicit none + private + public :: lammps_open, lammps_open_no_mpi, lammps_close, lammps_file, & + lammps_command, lammps_free, lammps_extract_global, & + lammps_extract_atom, lammps_extract_compute, lammps_extract_fix, & + lammps_extract_variable, lammps_get_natoms, lammps_gather_atoms, & + lammps_scatter_atoms + public :: lammps_instance, C_ptr, C_double, C_int + + !! Functions supplemental to the prototypes in library.h. {{{1 + !! The function definitions (in C++) are contained in LAMMPS-wrapper.cpp. + !! I would have written the first in Fortran, but the MPI libraries (which + !! were written in C) have C-based functions to convert from Fortran MPI + !! handles to C MPI handles, and there is no Fortran equivalent for those + !! functions. + interface + subroutine lammps_open_wrapper (argc, argv, communicator, ptr) & + bind (C, name='lammps_open_fortran_wrapper') + import :: C_int, C_ptr + integer (C_int), value :: argc + type (C_ptr), dimension(*) :: argv + integer, value :: communicator + type (C_ptr) :: ptr + end subroutine lammps_open_wrapper + subroutine lammps_actual_error_all (ptr, file, line, str) & + bind (C, name='lammps_error_all') + import :: C_int, C_char, C_ptr + type (C_ptr), value :: ptr + character (kind=C_char), dimension(*), intent(in) :: file, str + integer (C_int), value :: line + end subroutine lammps_actual_error_all + function lammps_get_ntypes (ptr) result (ntypes) & + bind (C, name='lammps_get_ntypes') + import :: C_int, C_ptr + type (C_ptr), value :: ptr + integer (C_int) :: ntypes + end function lammps_get_ntypes + function lammps_actual_extract_compute_vectorsize (ptr, id, style) & + result (vectorsize) bind (C, name='lammps_extract_compute_vectorsize') + import :: C_int, C_char, C_ptr + integer (C_int) :: vectorsize + type (C_ptr), value :: ptr + character (kind=C_char), dimension(*) :: id + integer (C_int), value :: style + end function lammps_actual_extract_compute_vectorsize + subroutine lammps_actual_extract_compute_arraysize (ptr, id, style, & + nrows, ncols) bind (C, name='lammps_extract_compute_arraysize') + import :: C_int, C_char, C_ptr + integer (C_int) :: arraysize + type (C_ptr), value :: ptr + character (kind=C_char), dimension(*) :: id + integer (C_int), value :: style + integer (C_int) :: nrows, ncols + end subroutine lammps_actual_extract_compute_arraysize + function lammps_actual_extract_fix_vectorsize (ptr, id, style) & + result (vectorsize) bind (C, name='lammps_extract_fix_vectorsize') + import :: C_int, C_char, C_ptr + integer (C_int) :: vectorsize + type (C_ptr), value :: ptr + character (kind=C_char), dimension(*) :: id + integer (C_int), value :: style + end function lammps_actual_extract_fix_vectorsize + subroutine lammps_actual_extract_fix_arraysize (ptr, id, style, & + nrows, ncols) bind (C, name='lammps_extract_fix_arraysize') + import :: C_int, C_char, C_ptr + type (C_ptr), value :: ptr + character (kind=C_char), dimension(*) :: id + integer (C_int), value :: style + integer (C_int) :: nrows, ncols + end subroutine lammps_actual_extract_fix_arraysize + end interface + + !! Functions/subroutines defined in library.h and library.cpp {{{1 + interface + subroutine lammps_actual_open_no_mpi (argc, argv, ptr) & + bind (C, name='lammps_open_no_mpi') + import :: C_int, C_ptr + integer (C_int), value :: argc + type (C_ptr), dimension(*) :: argv + type (C_ptr) :: ptr + end subroutine lammps_actual_open_no_mpi + + subroutine lammps_close (ptr) bind (C, name='lammps_close') + import :: C_ptr + type (C_ptr), value :: ptr + end subroutine lammps_close + + subroutine lammps_actual_file (ptr, str) bind (C, name='lammps_file') + import :: C_ptr, C_char + type (C_ptr), value :: ptr + character (kind=C_char), dimension(*) :: str + end subroutine lammps_actual_file + + function lammps_actual_command (ptr, str) result (command) & + bind (C, name='lammps_command') + import :: C_ptr, C_char + type (C_ptr), value :: ptr + character (kind=C_char), dimension(*) :: str + type (C_ptr) :: command + end function lammps_actual_command + + subroutine lammps_free (ptr) bind (C, name='lammps_free') + import :: C_ptr + type (C_ptr), value :: ptr + end subroutine lammps_free + + function lammps_actual_extract_global (ptr, name) & + bind (C, name='lammps_extract_global') result (global) + import :: C_ptr, C_char + type (C_ptr), value :: ptr + character (kind=C_char), dimension(*) :: name + type (C_ptr) :: global + end function lammps_actual_extract_global + + function lammps_actual_extract_atom (ptr, name) & + bind (C, name='lammps_extract_atom') result (atom) + import :: C_ptr, C_char + type (C_ptr), value :: ptr + character (kind=C_char), dimension(*) :: name + type (C_ptr) :: atom + end function lammps_actual_extract_atom + + function lammps_actual_extract_compute (ptr, id, style, type) & + result (compute) bind (C, name='lammps_extract_compute') + import :: C_ptr, C_char, C_int + type (C_ptr), value :: ptr + character (kind=C_char), dimension(*) :: id + integer (C_int), value :: style, type + type (C_ptr) :: compute + end function lammps_actual_extract_compute + + function lammps_actual_extract_fix (ptr, id, style, type, i, j) & + result (fix) bind (C, name='lammps_extract_fix') + import :: C_ptr, C_char, C_int + type (C_ptr), value :: ptr + character (kind=C_char), dimension(*) :: id + integer (C_int), value :: style, type, i, j + type (C_ptr) :: fix + end function lammps_actual_extract_fix + + function lammps_actual_extract_variable (ptr, name, group) & + result (variable) bind (C, name='lammps_extract_variable') + import :: C_ptr, C_char + type (C_ptr), value :: ptr + character (kind=C_char), dimension(*) :: name, group + type (C_ptr) :: variable + end function lammps_actual_extract_variable + + function lammps_get_natoms (ptr) result (natoms) & + bind (C, name='lammps_get_natoms') + import :: C_ptr, C_int + type (C_ptr), value :: ptr + integer (C_int) :: natoms + end function lammps_get_natoms + + subroutine lammps_actual_gather_atoms (ptr, name, type, count, data) & + bind (C, name='lammps_gather_atoms') + import :: C_ptr, C_int, C_char + type (C_ptr), value :: ptr, data + character (kind=C_char), dimension(*) :: name + integer (C_int), value :: type, count + end subroutine lammps_actual_gather_atoms + + subroutine lammps_actual_scatter_atoms (ptr, name, type, count, data) & + bind (C, name='lammps_scatter_atoms') + import :: C_ptr, C_int, C_char + type (C_ptr), value :: ptr, data + character (kind=C_char), dimension(*) :: name + integer (C_int), value :: type, count + end subroutine lammps_actual_scatter_atoms + end interface + + ! Generic functions for the wrappers below {{{1 + + interface lammps_extract_global + module procedure lammps_extract_global_i, & + lammps_extract_global_dp + end interface lammps_extract_global + + interface lammps_extract_atom + module procedure lammps_extract_atom_ia, & + lammps_extract_atom_dpa, & + lammps_extract_atom_dp2a + end interface lammps_extract_atom + + interface lammps_extract_compute + module procedure lammps_extract_compute_dp, & + lammps_extract_compute_dpa, & + lammps_extract_compute_dp2a + end interface lammps_extract_compute + + interface lammps_extract_fix + module procedure lammps_extract_fix_dp, & + lammps_extract_fix_dpa, & + lammps_extract_fix_dp2a + end interface lammps_extract_fix + + interface lammps_extract_variable + module procedure lammps_extract_variable_dp, & + lammps_extract_variable_dpa + end interface lammps_extract_variable + + interface lammps_gather_atoms + module procedure lammps_gather_atoms_ia, lammps_gather_atoms_dpa + end interface lammps_gather_atoms + + interface lammps_scatter_atoms + module procedure lammps_scatter_atoms_ia, lammps_scatter_atoms_dpa + end interface lammps_scatter_atoms + +contains !! Wrapper functions local to this module {{{1 + + subroutine lammps_open (command_line, communicator, ptr) + character (len=*), intent(in) :: command_line + integer, intent(in) :: communicator + type (C_ptr) :: ptr + integer (C_int) :: argc + type (C_ptr), dimension(:), allocatable :: argv + character (kind=C_char), dimension(len_trim(command_line)+1), target :: & + c_command_line + c_command_line = string2Cstring (command_line) + call Cstring2argcargv (c_command_line, argc, argv) + call lammps_open_wrapper (argc, argv, communicator, ptr) + deallocate (argv) + end subroutine lammps_open + +!----------------------------------------------------------------------------- + + subroutine lammps_open_no_mpi (command_line, ptr) + character (len=*), intent(in) :: command_line + type (C_ptr) :: ptr + integer (C_int) :: argc + type (C_ptr), dimension(:), allocatable :: argv + character (kind=C_char), dimension(len_trim(command_line)+1), target :: & + c_command_line + c_command_line = string2Cstring (command_line) + call Cstring2argcargv (c_command_line, argc, argv) + call lammps_actual_open_no_mpi (argc, argv, ptr) + deallocate (argv) + end subroutine lammps_open_no_mpi + +!----------------------------------------------------------------------------- + + subroutine lammps_file (ptr, str) + type (C_ptr) :: ptr + character (len=*) :: str + character (kind=C_char), dimension(len_trim(str)+1) :: Cstr + Cstr = string2Cstring (str) + call lammps_actual_file (ptr, Cstr) + end subroutine lammps_file + +!----------------------------------------------------------------------------- + + subroutine lammps_command (ptr, str) + type (C_ptr) :: ptr + character (len=*) :: str + character (kind=C_char), dimension(len_trim(str)+1) :: Cstr + type (C_ptr) :: dummy + Cstr = string2Cstring (str) + dummy = lammps_actual_command (ptr, Cstr) + end subroutine lammps_command + +!----------------------------------------------------------------------------- + +! lammps_extract_global {{{2 + function lammps_extract_global_Cptr (ptr, name) result (global) + type (C_ptr) :: global + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: name + character (kind=C_char), dimension(len_trim(name)+1) :: Cname + Cname = string2Cstring (name) + global = lammps_actual_extract_global (ptr, Cname) + end function lammps_extract_global_Cptr + subroutine lammps_extract_global_i (global, ptr, name) + integer (C_int), pointer, intent(out) :: global + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: name + type (C_ptr) :: Cptr + Cptr = lammps_extract_global_Cptr (ptr, name) + call C_F_pointer (Cptr, global) + end subroutine lammps_extract_global_i + subroutine lammps_extract_global_dp (global, ptr, name) + real (C_double), pointer, intent(out) :: global + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: name + type (C_ptr) :: Cptr + Cptr = lammps_extract_global_Cptr (ptr, name) + call C_F_pointer (Cptr, global) + end subroutine lammps_extract_global_dp + +!----------------------------------------------------------------------------- + +! lammps_extract_atom {{{2 + function lammps_extract_atom_Cptr (ptr, name) result (atom) + type (C_ptr) :: atom + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: name + character (kind=C_char), dimension(len_trim(name)+1) :: Cname + Cname = string2Cstring (name) + atom = lammps_actual_extract_atom (ptr, Cname) + end function lammps_extract_atom_Cptr + subroutine lammps_extract_atom_ia (atom, ptr, name) + integer (C_int), dimension(:), pointer, intent(out) :: atom + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: name + type (C_ptr) :: Cptr + integer (C_int), pointer :: nelements + call lammps_extract_global_i (nelements, ptr, 'nlocal') + Cptr = lammps_extract_atom_Cptr (ptr, name) + call C_F_pointer (Cptr, atom, (/nelements/)) + end subroutine lammps_extract_atom_ia + subroutine lammps_extract_atom_dpa (atom, ptr, name) + real (C_double), dimension(:), pointer, intent(out) :: atom + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: name + type (C_ptr) :: Cptr + integer (C_int), pointer :: nlocal + integer :: nelements + real (C_double), dimension(:), pointer :: Fptr + if ( name == 'mass' ) then + nelements = lammps_get_ntypes (ptr) + 1 + else if ( name == 'x' .or. name == 'v' .or. name == 'f' .or. & + name == 'mu' .or. name == 'omega' .or. name == 'torque' .or. & + name == 'angmom' ) then + ! We should not be getting a rank-2 array here! + call lammps_error_all (ptr, FLERR, 'You cannot extract those atom& + & data (' // trim(name) // ') into a rank 1 array.') + return + else + ! Everything else we can get is probably nlocal units long + call lammps_extract_global_i (nlocal, ptr, 'nlocal') + nelements = nlocal + end if + Cptr = lammps_extract_atom_Cptr (ptr, name) + call C_F_pointer (Cptr, Fptr, (/nelements/)) + if ( name == 'mass' ) then + atom(0:) => Fptr + else + atom => Fptr + end if + end subroutine lammps_extract_atom_dpa + subroutine lammps_extract_atom_dp2a (atom, ptr, name) + real (C_double), dimension(:,:), pointer, intent(out) :: atom + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: name + type (C_ptr) :: Cptr + type (C_ptr), pointer, dimension(:) :: Catom + integer (C_int), pointer :: nelements + if ( name /= 'x' .and. name /= 'v' .and. name /= 'f' .and. & + name /= 'mu' .and. name /= 'omega' .and. name /= 'tandque' .and. & + name /= 'angmom' ) then + ! We should not be getting a rank-2 array here! + call lammps_error_all (ptr, FLERR, 'You cannot extract those atom& + & data (' // trim(name) // ') into a rank 2 array.') + return + end if + Cptr = lammps_extract_atom_Cptr (ptr, name) + call lammps_extract_global_i (nelements, ptr, 'nlocal') + ! Catom will now be the array of void* pointers that the void** pointer + ! pointed to. Catom(1) is now the pointer to the first element. + call C_F_pointer (Cptr, Catom, (/nelements/)) + ! Now get the actual array, which has its shape transposed from what we + ! might think of it in C + call C_F_pointer (Catom(1), atom, (/3, nelements/)) + end subroutine lammps_extract_atom_dp2a + +!----------------------------------------------------------------------------- + +! lammps_extract_compute {{{2 + function lammps_extract_compute_Cptr (ptr, id, style, type) result (compute) + type (C_ptr) :: compute + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: id + integer, intent(in) :: style, type + integer (kind=C_int) :: Cstyle, Ctype + character (kind=C_char), dimension(len_trim(id)+1) :: Cid + Cid = string2Cstring (id) + Cstyle = style + Ctype = type + compute = lammps_actual_extract_compute (ptr, Cid, Cstyle, Ctype) + end function lammps_extract_compute_Cptr + subroutine lammps_extract_compute_dp (compute, ptr, id, style, type) + real (C_double), pointer, intent(out) :: compute + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: id + integer, intent(in) :: style, type + type (C_ptr) :: Cptr + ! The only valid values of (style,type) are (0,0) for scalar 'compute' + if ( style /= 0 ) then + call lammps_error_all (ptr, FLERR, 'You cannot pack per-atom/local& + & data into a scalar.') + return + end if + if ( type == 1 ) then + call lammps_error_all (ptr, FLERR, 'You cannot extract a compute& + & vector (rank 1) into a scalar.') + return + else if ( type == 2 ) then + call lammps_error_all (ptr, FLERR, 'You cannot extract a compute& + & array (rank 2) into a scalar.') + return + end if + Cptr = lammps_extract_compute_Cptr (ptr, id, style, type) + call C_F_pointer (Cptr, compute) + end subroutine lammps_extract_compute_dp + subroutine lammps_extract_compute_dpa (compute, ptr, id, style, type) + real (C_double), dimension(:), pointer, intent(out) :: compute + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: id + integer, intent(in) :: style, type + type (C_ptr) :: Cptr + integer :: nelements + ! Check for the correct dimensionality + if ( type == 0 ) then + call lammps_error_all (ptr, FLERR, 'You cannot extract a compute& + & scalar (rank 0) into a rank 1 variable.') + return + else if ( type == 2 ) then + call lammps_error_all (ptr, FLERR, 'You cannot extract a compute& + & array (rank 2) into a rank 1 variable.') + return + end if + nelements = lammps_extract_compute_vectorsize (ptr, id, style) + Cptr = lammps_extract_compute_Cptr (ptr, id, style, type) + call C_F_pointer (Cptr, compute, (/nelements/)) + end subroutine lammps_extract_compute_dpa + subroutine lammps_extract_compute_dp2a (compute, ptr, id, style, type) + real (C_double), dimension(:,:), pointer, intent(out) :: compute + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: id + integer, intent(in) :: style, type + type (C_ptr) :: Cptr + type (C_ptr), pointer, dimension(:) :: Ccompute + integer :: nr, nc + ! Check for the correct dimensionality + if ( type == 0 ) then + call lammps_error_all (ptr, FLERR, 'You cannot extract a compute& + & scalar (rank 0) into a rank 2 variable.') + return + else if ( type == 1 ) then + call lammps_error_all (ptr, FLERR, 'You cannot extract a compute& + & array (rank 1) into a rank 2 variable.') + return + end if + call lammps_extract_compute_arraysize (ptr, id, style, nr, nc) + call C_F_pointer (Cptr, Ccompute, (/nr/)) + ! Note that the matrix is transposed, from Fortran's perspective + call C_F_pointer (Ccompute(1), compute, (/nc, nr/)) + end subroutine lammps_extract_compute_dp2a + +!----------------------------------------------------------------------------- + +! lammps_extract_fix {{{2 + function lammps_extract_fix_Cptr (ptr, id, style, type, i, j) & + result (fix) + type (C_ptr) :: fix + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: id + integer, intent(in) :: style, type, i, j + character (kind=C_char), dimension(len_trim(id)+1) :: Cid + integer (kind=C_int) :: Cstyle, Ctype, Ci, Cj + Cid = string2Cstring (id) + Cstyle = style + Ctype = type + Ci = i - 1 ! This is for consistency with the values from f_ID[i], + Cj = j - 1 ! which is different from what library.cpp uses! + if ( (type >= 1 .and. Ci < 0) .or. & + (type == 2 .and. (Ci < 0 .or. Cj < 0) ) ) then + call lammps_error_all (ptr, FLERR, 'Index out of range in& + & lammps_extract_fix') + end if + fix = lammps_actual_extract_fix (ptr, Cid, Cstyle, Ctype, Ci, Cj) + end function lammps_extract_fix_Cptr + subroutine lammps_extract_fix_dp (fix, ptr, id, style, type, i, j) + real (C_double), intent(out) :: fix + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: id + integer, intent(in) :: style, type, i, j + type (C_ptr) :: Cptr + real (C_double), pointer :: Fptr + ! Check for the correct dimensionality + if ( style /= 0 ) then + select case (type) + case (0) + call lammps_error_all (ptr, FLERR, 'There is no per-atom or local& + & scalar data available from fixes.') + case (1) + call lammps_error_all (ptr, FLERR, 'You cannot extract a fix''s & + &per-atom/local vector (rank 1) into a scalar.') + case (2) + call lammps_error_all (ptr, FLERR, 'You cannot extract a fix''s & + &per-atom/local array (rank 2) into a scalar.') + case default + call lammps_error_all (ptr, FLERR, 'Invalid extract_fix style/& + &type combination.') + end select + return + end if + Cptr = lammps_extract_fix_Cptr (ptr, id, style, type, i, j) + call C_F_pointer (Cptr, Fptr) + fix = Fptr + nullify (Fptr) + ! Memory is only allocated for "global" fix variables + if ( style == 0 ) call lammps_free (Cptr) + end subroutine lammps_extract_fix_dp + subroutine lammps_extract_fix_dpa (fix, ptr, id, style, type, i, j) + real (C_double), dimension(:), pointer, intent(out) :: fix + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: id + integer, intent(in) :: style, type, i, j + type (C_ptr) :: Cptr + integer :: fix_len + ! Check for the correct dimensionality + if ( style == 0 ) then + call lammps_error_all (ptr, FLERR, 'You can''t extract the& + & whole vector from global fix data') + return + else if ( type == 0 ) then + call lammps_error_all (ptr, FLERR, 'You can''t extract a fix& + & scalar into a rank 1 variable') + return + else if ( type == 2 ) then + call lammps_error_all (ptr, FLERR, 'You cannot extract a fix& + & array into a rank 1 variable.') + return + else if ( type /= 1 ) then + call lammps_error_all (ptr, FLERR, 'Invalid type for fix extraction.') + return + end if + fix_len = lammps_extract_fix_vectorsize (ptr, id, style) + call C_F_pointer (Cptr, fix, (/fix_len/)) + ! Memory is only allocated for "global" fix variables, which we should + ! never get here, so no need to call lammps_free! + end subroutine lammps_extract_fix_dpa + subroutine lammps_extract_fix_dp2a (fix, ptr, id, style, type, i, j) + real (C_double), dimension(:,:), pointer, intent(out) :: fix + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: id + integer, intent(in) :: style, type, i, j + type (C_ptr) :: Cptr + type (C_ptr), pointer, dimension(:) :: Cfix + integer :: nr, nc + ! Check for the correct dimensionality + if ( style == 0 ) then + call lammps_error_all (ptr, FLERR, 'It is not possible to extract the& + & entire array from global fix data.') + return + else if ( type == 0 ) then + call lammps_error_all (ptr, FLERR, 'You cannot extract a fix& + & scalar (rank 0) into a rank 2 variable.') + return + else if ( type == 1 ) then + call lammps_error_all (ptr, FLERR, 'You cannot extract a fix& + & vector (rank 1) into a rank 2 variable.') + return + end if + call lammps_extract_fix_arraysize (ptr, id, style, nr, nc) + ! Extract pointer to first element as Cfix(1) + call C_F_pointer (Cptr, Cfix, (/nr/)) + ! Now extract the array, which is transposed + call C_F_pointer (Cfix(1), fix, (/nc, nr/)) + end subroutine lammps_extract_fix_dp2a + +!----------------------------------------------------------------------------- + +! lammps_extract_variable {{{2 + function lammps_extract_variable_Cptr (ptr, name, group) result (variable) + type (C_ptr) :: ptr, variable + character (len=*) :: name + character (len=*), optional :: group + character (kind=C_char), dimension(len_trim(name)+1) :: Cname + character (kind=C_char), dimension(:), allocatable :: Cgroup + Cname = string2Cstring (name) + if ( present(group) ) then + allocate (Cgroup(len_trim(group)+1)) + Cgroup = string2Cstring (group) + else + allocate (Cgroup(1)) + Cgroup(1) = C_NULL_CHAR + end if + variable = lammps_actual_extract_variable (ptr, Cname, Cgroup) + deallocate (Cgroup) + end function lammps_extract_variable_Cptr + subroutine lammps_extract_variable_dp (variable, ptr, name, group) + real (C_double), intent(out) :: variable + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: name + character (len=*), intent(in), optional :: group + type (C_ptr) :: Cptr + real (C_double), pointer :: Fptr + if ( present(group) ) then + Cptr = lammps_extract_variable_Cptr (ptr, name, group) + else + Cptr = lammps_extract_variable_Cptr (ptr, name) + end if + call C_F_pointer (Cptr, Fptr) + variable = Fptr + nullify (Fptr) + call lammps_free (Cptr) + end subroutine lammps_extract_variable_dp + subroutine lammps_extract_variable_dpa (variable, ptr, name, group) + real (C_double), dimension(:), allocatable, intent(out) :: variable + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: name + character (len=*), intent(in), optional :: group + type (C_ptr) :: Cptr + real (C_double), dimension(:), pointer :: Fptr + integer :: natoms + if ( present(group) ) then + Cptr = lammps_extract_variable_Cptr (ptr, name, group) + else + Cptr = lammps_extract_variable_Cptr (ptr, name) + end if + natoms = lammps_get_natoms (ptr) + allocate (variable(natoms)) + call C_F_pointer (Cptr, Fptr, (/natoms/)) + variable = Fptr + nullify (Fptr) + call lammps_free (Cptr) + end subroutine lammps_extract_variable_dpa + +!-------------------------------------------------------------------------2}}} + + subroutine lammps_gather_atoms_ia (ptr, name, count, data) + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: name + integer, intent(in) :: count + integer, dimension(:), allocatable, intent(out) :: data + type (C_ptr) :: Cdata + integer (C_int), dimension(:), pointer :: Fdata + integer (C_int) :: natoms + character (kind=C_char), dimension(len_trim(name)+1) :: Cname + integer (C_int), parameter :: Ctype = 0_C_int + integer (C_int) :: Ccount + natoms = lammps_get_natoms (ptr) + Cname = string2Cstring (name) + if ( count /= 1 .and. count /= 3 ) then + call lammps_error_all (ptr, FLERR, 'lammps_gather_atoms requires& + & count to be either 1 or 3') + else + Ccount = count + end if + allocate ( Fdata(count*natoms) ) + allocate ( data(count*natoms) ) + Cdata = C_loc (Fdata(1)) + call lammps_actual_gather_atoms (ptr, Cname, Ctype, Ccount, Cdata) + data = Fdata + deallocate (Fdata) + end subroutine lammps_gather_atoms_ia + subroutine lammps_gather_atoms_dpa (ptr, name, count, data) + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: name + integer, intent(in) :: count + double precision, dimension(:), allocatable, intent(out) :: data + type (C_ptr) :: Cdata + real (C_double), dimension(:), pointer :: Fdata + integer (C_int) :: natoms + character (kind=C_char), dimension(len_trim(name)+1) :: Cname + integer (C_int), parameter :: Ctype = 1_C_int + integer (C_int) :: Ccount + natoms = lammps_get_natoms (ptr) + Cname = string2Cstring (name) + if ( count /= 1 .and. count /= 3 ) then + call lammps_error_all (ptr, FLERR, 'lammps_gather_atoms requires& + & count to be either 1 or 3') + else + Ccount = count + end if + allocate ( Fdata(count*natoms) ) + allocate ( data(count*natoms) ) + Cdata = C_loc (Fdata(1)) + call lammps_actual_gather_atoms (ptr, Cname, Ctype, Ccount, Cdata) + data = Fdata(:) + deallocate (Fdata) + end subroutine lammps_gather_atoms_dpa + +!----------------------------------------------------------------------------- + + subroutine lammps_scatter_atoms_ia (ptr, name, data) + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: name + integer, dimension(:), intent(in) :: data + integer (kind=C_int) :: natoms, Ccount + integer (kind=C_int), parameter :: Ctype = 0_C_int + character (kind=C_char), dimension(len_trim(name)+1) :: Cname + integer (C_int), dimension(size(data)), target :: Fdata + type (C_ptr) :: Cdata + natoms = lammps_get_natoms (ptr) + Cname = string2Cstring (name) + Ccount = size(data) / natoms + if ( Ccount /= 1 .and. Ccount /= 3 ) & + call lammps_error_all (ptr, FLERR, 'lammps_gather_atoms requires& + & count to be either 1 or 3') + Fdata = data + Cdata = C_loc (Fdata(1)) + call lammps_actual_scatter_atoms (ptr, Cname, Ctype, Ccount, Cdata) + end subroutine lammps_scatter_atoms_ia + subroutine lammps_scatter_atoms_dpa (ptr, name, data) + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: name + double precision, dimension(:), intent(in) :: data + integer (kind=C_int) :: natoms, Ccount + integer (kind=C_int), parameter :: Ctype = 1_C_int + character (kind=C_char), dimension(len_trim(name)+1) :: Cname + real (C_double), dimension(size(data)), target :: Fdata + type (C_ptr) :: Cdata + natoms = lammps_get_natoms (ptr) + Cname = string2Cstring (name) + Ccount = size(data) / natoms + if ( Ccount /= 1 .and. Ccount /= 3 ) & + call lammps_error_all (ptr, FLERR, 'lammps_gather_atoms requires& + & count to be either 1 or 3') + Fdata = data + Cdata = C_loc (Fdata(1)) + call lammps_actual_scatter_atoms (ptr, Cname, Ctype, Ccount, Cdata) + end subroutine lammps_scatter_atoms_dpa + +!----------------------------------------------------------------------------- + + function lammps_extract_compute_vectorsize (ptr, id, style) & + result (vectorsize) + integer :: vectorsize + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: id + integer, intent(in) :: style + integer (C_int) :: Cvectorsize, Cstyle + character (kind=C_char), dimension(len_trim(id)+1) :: Cid + Cid = string2Cstring (id) + Cstyle = int(style, C_int) + Cvectorsize = lammps_actual_extract_compute_vectorsize (ptr, Cid, Cstyle) + vectorsize = int(Cvectorsize, kind(vectorsize)) + end function lammps_extract_compute_vectorsize + +!----------------------------------------------------------------------------- + + function lammps_extract_fix_vectorsize (ptr, id, style) & + result (vectorsize) + integer :: vectorsize + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: id + integer, intent(in) :: style + integer (C_int) :: Cvectorsize, Cstyle + character (kind=C_char), dimension(len_trim(id)+1) :: Cid + Cid = string2Cstring (id) + Cstyle = int(style, C_int) + Cvectorsize = lammps_actual_extract_fix_vectorsize (ptr, Cid, Cstyle) + vectorsize = int(Cvectorsize, kind(vectorsize)) + end function lammps_extract_fix_vectorsize + +!----------------------------------------------------------------------------- + + subroutine lammps_extract_compute_arraysize (ptr, id, style, nrows, ncols) + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: id + integer, intent(in) :: style + integer, intent(out) :: nrows, ncols + integer (C_int) :: Cstyle, Cnrows, Cncols + character (kind=C_char), dimension(len_trim(id)+1) :: Cid + Cid = string2Cstring (id) + Cstyle = int (style, C_int) + call lammps_actual_extract_compute_arraysize (ptr, Cid, Cstyle, & + Cnrows, Cncols) + nrows = int (Cnrows, kind(nrows)) + ncols = int (Cncols, kind(ncols)) + end subroutine lammps_extract_compute_arraysize + +!----------------------------------------------------------------------------- + + subroutine lammps_extract_fix_arraysize (ptr, id, style, nrows, ncols) + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: id + integer, intent(in) :: style + integer, intent(out) :: nrows, ncols + integer (C_int) :: Cstyle, Cnrows, Cncols + character (kind=C_char), dimension(len_trim(id)+1) :: Cid + Cid = string2Cstring (id) + Cstyle = int (style, kind(Cstyle)) + call lammps_actual_extract_fix_arraysize (ptr, Cid, Cstyle, & + Cnrows, Cncols) + nrows = int (Cnrows, kind(nrows)) + ncols = int (Cncols, kind(ncols)) + end subroutine lammps_extract_fix_arraysize + +!----------------------------------------------------------------------------- + + subroutine lammps_error_all (ptr, file, line, str) + type (C_ptr), intent(in) :: ptr + character (len=*), intent(in) :: file, str + integer, intent(in) :: line + character (kind=C_char), dimension(len_trim(file)+1) :: Cfile + character (kind=C_char), dimension(len_trim(str)+1) :: Cstr + integer (C_int) :: Cline + Cline = int(line, kind(Cline)) + Cfile = string2Cstring (file) + Cstr = string2Cstring (str) + call lammps_actual_error_all (ptr, Cfile, Cline, Cstr) + end subroutine lammps_error_all + +!----------------------------------------------------------------------------- + +! Locally defined helper functions {{{1 + + pure function string2Cstring (string) result (C_string) + use, intrinsic :: ISO_C_binding, only : C_char, C_NULL_CHAR + character (len=*), intent(in) :: string + character (len=1, kind=C_char) :: C_string (len_trim(string)+1) + integer :: i, n + n = len_trim (string) + forall (i = 1:n) + C_string(i) = string(i:i) + end forall + C_string(n+1) = C_NULL_CHAR + end function string2Cstring + +!----------------------------------------------------------------------------- + + subroutine Cstring2argcargv (Cstring, argc, argv) + !! Converts a C-style string to argc and argv, that is, words in Cstring + !! become C-style strings in argv. IMPORTANT: Cstring is modified by + !! this routine! I would make Cstring local TO this routine and accept + !! a Fortran-style string instead, but we run into scoping and + !! allocation problems that way. This routine assumes the string is + !! null-terminated, as all C-style strings must be. + + character (kind=C_char), dimension(*), target, intent(inout) :: Cstring + integer (C_int), intent(out) :: argc + type (C_ptr), dimension(:), allocatable, intent(out) :: argv + + integer :: StringStart, SpaceIndex, strlen, argnum + + argc = 1_C_int + + ! Find the length of the string + strlen = 1 + do while ( Cstring(strlen) /= C_NULL_CHAR ) + strlen = strlen + 1 + end do + + ! Find the number of non-escaped spaces + SpaceIndex = 2 + do while ( SpaceIndex < strlen ) + if ( Cstring(SpaceIndex) == ' ' .and. & + Cstring(SpaceIndex-1) /= '\' ) then + argc = argc + 1_C_int + ! Find the next non-space character + do while ( Cstring(SpaceIndex+1) == ' ') + SpaceIndex = SpaceIndex + 1 + end do + end if + SpaceIndex = SpaceIndex + 1 + end do + + ! Now allocate memory for argv + allocate (argv(argc)) + + ! Now find the string starting and ending locations + StringStart = 1 + SpaceIndex = 2 + argnum = 1 + do while ( SpaceIndex < strlen ) + if ( Cstring(SpaceIndex) == ' ' .and. & + Cstring(SpaceIndex-1) /= '\' ) then + ! Found a real space => split strings and store this one + Cstring(Spaceindex) = C_NULL_CHAR ! Replaces space with NULL + argv(argnum) = C_loc(Cstring(StringStart)) + argnum = argnum + 1 + ! Find the next non-space character + do while ( Cstring(SpaceIndex+1) == ' ') + SpaceIndex = SpaceIndex + 1 + end do + StringStart = SpaceIndex + 1 + else if ( Cstring(SpaceIndex) == ' ' .and. & + Cstring(SpaceIndex-1) == '\' ) then + ! Escaped space => remove backslash and move rest of array + Cstring(SpaceIndex-1:strlen-1) = Cstring(SpaceIndex:strlen) + strlen = strlen - 1 ! Last character is still C_NULL_CHAR + end if + SpaceIndex = SpaceIndex + 1 + end do + ! Now handle the last argument + argv(argnum) = C_loc(Cstring(StringStart)) + + end subroutine Cstring2argcargv + +! 1}}} + +end module LAMMPS + +! vim: foldmethod=marker tabstop=3 softtabstop=3 shiftwidth=3 expandtab diff --git a/examples/COUPLE/fortran2/README b/examples/COUPLE/fortran2/README index 03c76fb3ca..01eb76b0a1 100644 --- a/examples/COUPLE/fortran2/README +++ b/examples/COUPLE/fortran2/README @@ -1,221 +1,265 @@ -LAMMPS.F90 defines a Fortran 2003 module, LAMMPS, which wraps all functions in -src/library.h so they can be used directly from Fortran-encoded programs. - -All functions in src/library.h that use and/or return C-style pointers have -Fortran wrapper functions that use Fortran-style arrays, pointers, and -strings; all C-style memory management is handled internally with no user -intervention. - -This interface was created by Karl Hammond who you can contact with -questions: - -Karl D. Hammond -University of Tennessee, Knoxville -karlh at ugcs.caltech.edu -karlh at utk.edu - -------------------------------------- - ---COMPILATION-- - -First, be advised that mixed-language programming is not trivial. It requires -you to link in the required libraries of all languages you use (in this case, -those for Fortran, C, and C++), as well as any other libraries required. -You are also advised to read the --USE-- section below before trying to -compile. - -The following steps will work to compile this module (replace ${LAMMPS_SRC} -with the path to your LAMMPS source directory): - (1) Compile LAMMPS as a static library. Call the resulting file ${LAMMPS_LIB}, - which will have an actual name lake liblmp_openmpi.a. If compiling - using the MPI stubs in ${LAMMPS_SRC}/STUBS, you will need to know where - libmpi.a is as well (I'll call it ${MPI_STUBS} hereafter) - (2) Copy said library to your Fortran program's source directory or include - its location in a -L${LAMMPS_SRC} flag to your compiler. - (3) Compile (but don't link!) LAMMPS.F90. Example: - mpif90 -c LAMMPS.f90 - OR - gfortran -c LAMMPS.F90 - Copy the LAMMPS.o and lammps.mod (or whatever your compiler calls module - files) to your Fortran program's source directory. - NOTE: you may get a warning such as, - subroutine lammps_open_wrapper (argc, argv, communicator, ptr) & - Variable 'communicator' at (1) is a parameter to the BIND(C) - procedure 'lammps_open_wrapper' but may not be C interoperable - This is normal (see --IMPLEMENTATION NOTES--). - (4) Compile (but don't link) LAMMPS-wrapper.cpp. You will need its header - file as well. You will have to provide the locations of LAMMPS's - header files. For example, - mpicxx -c -I${LAMMPS_SRC} LAMMPS-wrapper.cpp - OR - g++ -c -I${LAMMPS_SRC} -I${LAMMPS_SRC}/STUBS LAMMPS-wrapper.cpp - OR - icpc -c -I${LAMMPS_SRC} -I${LAMMPS_SRC}/STUBS LAMMPS-wrapper.cpp - Copy the resulting object file LAMMPS-wrapper.o to your Fortran program's - source directory. - (4b) OPTIONAL: Make a library so you can carry around two files instead of - three. Example: - ar rs liblammps_fortran.a LAMMPS.o LAMMPS-wrapper.o - This will create the file liblammps_fortran.a that you can use in place - of "LAMMPS.o LAMMPS-wrapper.o" in part (6). Note that you will still - need to have the .mod file from part (3). - - It is also possible to add LAMMPS.o and LAMMPS-wrapper.o into the - LAMMPS library (e.g., liblmp_openmpi.a) instead of creating a separate - library, like so: - ar rs ${LAMMPS_LIB} LAMMPS.o LAMMPS-wrapper.o - In this case, you can now use the Fortran wrapper functions as if they - were part of the usual LAMMPS library interface (if you have the module - file visible to the compiler, that is). - (5) Compile your Fortran program. Example: - mpif90 -c myfreeformatfile.f90 - mpif90 -c myfixedformatfile.f - OR - gfortran -c myfreeformatfile.f90 - gfortran -c myfixedformatfile.f - The object files generated by these steps are collectively referred to - as ${my_object_files} in the next step(s). - - IMPORTANT: If the Fortran module from part (3) is not in the current - directory or in one searched by the compiler for module files, you will - need to include that location via the -I flag to the compiler. - (6) Link everything together, including any libraries needed by LAMMPS (such - as the C++ standard library, the C math library, the JPEG library, fftw, - etc.) For example, - mpif90 LAMMPS.o LAMMPS-wrapper.o ${my_object_files} \ - ${LAMMPS_LIB} -lstdc++ -lm - OR - gfortran LAMMPS.o LAMMPS-wrapper.o ${my_object_files} \ - ${LAMMPS_LIB} ${MPI_STUBS} -lstdc++ -lm - OR - ifort LAMMPS.o LAMMPS-wrapper.o ${my_object_files} \ - ${LAMMPS_LIB} ${MPI_STUBS} -cxxlib -limf -lm - Any other required libraries (e.g. -ljpeg, -lfftw) should be added to - the end of this line. - -You should now have a working executable. - -Steps 3 and 4 above are accomplished, possibly after some modifications to -the makefile, by make using the attached makefile. - -------------------------------------- - ---USAGE-- - -To use this API, your program unit (PROGRAM/SUBROUTINE/FUNCTION/MODULE/etc.) -should look something like this: - program call_lammps - use LAMMPS - ! Other modules, etc. - implicit none - type (lammps_instance) :: lmp ! This is a pointer to your LAMMPS instance - double precision :: fix - double precision, dimension(:), allocatable :: fix2 - ! Rest of declarations - call lammps_open_no_mpi ('lmp -in /dev/null -screen out.lammps',lmp) - ! Set up rest of program here - call lammps_file (lmp, 'in.example') - call lammps_extract_fix (fix, lmp, '2', 0, 1, 1, 1) - call lammps_extract_fix (fix2, lmp, '4', 0, 2, 1, 1) - call lammps_close (lmp) - end program call_lammps - -Important notes: - * All arguments which are char* variables in library.cpp are character (len=*) - variables here. For example, - call lammps_command (lmp, 'units metal') - will work as expected. - * The public functions (the only ones you can use) have interfaces as - described in the comments at the top of LAMMPS.F90. They are not always - the same as those in library.h, since C strings are replaced by Fortran - strings and the like. - * The module attempts to check whether you have done something stupid (such - as assign a 2D array to a scalar), but it's not perfect. For example, the - command - call lammps_extract_global (nlocal, ptr, 'nlocal') - will give nlocal correctly if nlocal is of type INTEGER, but it will give - the wrong answer if nlocal is of type REAL or DOUBLE PRECISION. This is a - feature of the (void*) type cast in library.cpp. There is no way I can - check this for you! - * You are allowed to use REAL or DOUBLE PRECISION floating-point numbers. - All LAMMPS data (which are of type REAL(C_double)) are rounded off if - placed in single precision variables. It is tacitly assumed that NO C++ - variables are of type float; everything is int or double (since this is - all library.cpp currently handles). - * An example of a complete program is offered at the end of this file. - -------------------------------------- - ---TROUBLESHOOTING-- - -Compile-time errors probably indicate that your compiler is not new enough to -support Fortran 2003 features. For example, GCC 4.1.2 will not compile this -module, but GCC 4.4.0 will. - -If your compiler balks at 'use, intrinsic :: ISO_C_binding,' try removing the -intrinsic part so it looks like an ordinary module. However, it is likely -that such a compiler will also have problems with everything else in the -file as well. - -If you get a segfault as soon as the lammps_open call is made, check that you -compiled your program AND LAMMPS-header.cpp using the same MPI headers. Using -the stubs for one and the actual MPI library for the other will cause major -problems. - -If you find run-time errors, please pass them along via the LAMMPS Users -mailing list. Please provide a minimal working example along with the names -and versions of the compilers you are using. Please make sure the error is -repeatable and is in MY code, not yours (generating a minimal working example -will usually ensure this anyway). - -------------------------------------- - ---IMPLEMENTATION NOTES-- - -The Fortran procedures have the same names as the C procedures, and -their purpose is the same, but they may take different arguments. Here are -some of the important differences: - * lammps_open and lammps_open_no_mpi take a string instead of argc and - argv. This is necessary because C and C++ have a very different way - of treating strings than Fortran. - * All C++ functions that accept char* pointers now accept Fortran-style - strings within this interface instead. - * All of the lammps_extract_[something] functions, which return void* - C-style pointers, have been replaced by generic subroutines that return - Fortran variables (which may be arrays). The first argument houses the - variable to be returned; all other arguments are identical except as - stipulated above. Note that it is not possible to declare generic - functions that are selected based solely on the type/kind/rank (TKR) - signature of the return value, only based on the TKR of the arguments. - * The SHAPE of the first argument to lammps_extract_[something] is checked - against the "shape" of the C array (e.g., double vs. double* vs. double**). - Calling a subroutine with arguments of inappropriate rank will result in an - error at run time. - * All arrays passed to subroutines must be ALLOCATABLE and are REALLOCATED - to fit the shape of the array LAMMPS will be returning. - * The indices i and j in lammps_extract_fix are used the same way they - are in f_ID[i][j] references in LAMMPS (i.e., starting from 1). This is - different than the way library.cpp uses these numbers, but is more - consistent with the way arrays are accessed in LAMMPS and in Fortran. - * The char* pointer normally returned by lammps_command is thrown away - in this version; note also that lammps_command is now a subroutine - instead of a function. - * The pointer to LAMMPS itself is of type(lammps_instance), which is itself - a synonym for type(C_ptr), part of ISO_C_BINDING. Type (C_ptr) is - C's void* data type. This should be the only C data type that needs to - be used by the end user. - * This module will almost certainly generate a compile-time warning, - such as, - subroutine lammps_open_wrapper (argc, argv, communicator, ptr) & - Variable 'communicator' at (1) is a parameter to the BIND(C) - procedure 'lammps_open_wrapper' but may not be C interoperable - This happens because lammps_open_wrapper actually takes a Fortran - INTEGER argument, whose type is defined by the MPI library itself. The - Fortran integer is converted to a C integer by the MPI library (if such - conversion is actually necessary). - * Unlike library.cpp, this module returns COPIES of the data LAMMPS actually - uses. This is done for safety reasons, as you should, in general, not be - overwriting LAMMPS data directly from Fortran. If you require this - functionality, it is possible to write another function that, for example, - returns a Fortran pointer that resolves to the C/C++ data instead of - copying the contents of that pointer to the original array as is done now. +LAMMPS.F90 defines a Fortran 2003 module, LAMMPS, which wraps all functions in +src/library.h so they can be used directly from Fortran-encoded programs. + +All functions in src/library.h that use and/or return C-style pointers have +Fortran wrapper functions that use Fortran-style arrays, pointers, and +strings; all C-style memory management is handled internally with no user +intervention. See --USE-- for notes on how this interface differs from the +C interface (and the Python interface). + +This interface was created by Karl Hammond who you can contact with +questions: + +Karl D. Hammond +University of Tennessee, Knoxville +karlh at ugcs.caltech.edu +karlh at utk.edu + +------------------------------------- + +--COMPILATION-- + +First, be advised that mixed-language programming is not trivial. It requires +you to link in the required libraries of all languages you use (in this case, +those for Fortran, C, and C++), as well as any other libraries required. +You are also advised to read the --USE-- section below before trying to +compile. + +The following steps will work to compile this module (replace ${LAMMPS_SRC} +with the path to your LAMMPS source directory). + +Steps 3-5 are accomplished, possibly after some modifications to +the makefile, by make using the attached makefile. Said makefile also builds +the dynamically-linkable library (liblammps_fortran.so). + +** STATIC LIBRARY INSTRUCTIONS ** + (1) Compile LAMMPS as a static library. + Call the resulting file ${LAMMPS_LIB}, which will have an actual name + like liblmp_openmpi.a. If compiling using the MPI stubs in + ${LAMMPS_SRC}/STUBS, you will need to know where libmpi_stubs.a + is as well (I'll call it ${MPI_STUBS} hereafter) + (2) Copy said library to your Fortran program's source directory or replace + ${LAMMPS_LIB} with its full path in the instructions below. + (3) Compile (but don't link!) LAMMPS.F90. Example: + mpif90 -c LAMMPS.f90 + OR + gfortran -c LAMMPS.F90 + NOTE: you may get a warning such as, + subroutine lammps_open_wrapper (argc, argv, communicator, ptr) & + Variable 'communicator' at (1) is a parameter to the BIND(C) + procedure 'lammps_open_wrapper' but may not be C interoperable + This is normal (see --IMPLEMENTATION NOTES--). + + (4) Compile (but don't link) LAMMPS-wrapper.cpp. You will need its header + file as well. You will have to provide the locations of LAMMPS's + header files. For example, + mpicxx -c -I${LAMMPS_SRC} LAMMPS-wrapper.cpp + OR + g++ -c -I${LAMMPS_SRC} -I${LAMMPS_SRC}/STUBS LAMMPS-wrapper.cpp + OR + icpc -c -I${LAMMPS_SRC} -I${LAMMPS_SRC}/STUBS LAMMPS-wrapper.cpp + (5) OPTIONAL: Make a library from the object files so you can carry around + two files instead of three. Example: + ar rs liblammps_fortran.a LAMMPS.o LAMMPS-wrapper.o + This will create the file liblammps_fortran.a that you can use in place + of "LAMMPS.o LAMMPS-wrapper.o" later. Note that you will still + need to have the .mod file from part (3). + + It is also possible to add LAMMPS.o and LAMMPS-wrapper.o into the + LAMMPS library (e.g., liblmp_openmpi.a) instead of creating a separate + library, like so: + ar rs ${LAMMPS_LIB} LAMMPS.o LAMMPS-wrapper.o + In this case, you can now use the Fortran wrapper functions as if they + were part of the usual LAMMPS library interface (if you have the module + file visible to the compiler, that is). + (6) Compile (but don't link) your Fortran program. Example: + mpif90 -c myfreeformatfile.f90 + mpif90 -c myfixedformatfile.f + OR + gfortran -c myfreeformatfile.f90 + gfortran -c myfixedformatfile.f + The object files generated by these steps are collectively referred to + as ${my_object_files} in the next step(s). + + IMPORTANT: If the Fortran module from part (3) is not in the current + directory or in one searched by the compiler for module files, you will + need to include that location via the -I flag to the compiler, like so: + mpif90 -I${LAMMPS_SRC}/examples/COUPLE/fortran2 -c myfreeformatfile.f90 + + (7) Link everything together, including any libraries needed by LAMMPS (such + as the C++ standard library, the C math library, the JPEG library, fftw, + etc.) For example, + mpif90 LAMMPS.o LAMMPS-wrapper.o ${my_object_files} \ + ${LAMMPS_LIB} -lmpi_cxx -lstdc++ -lm + OR + gfortran LAMMPS.o LAMMPS-wrapper.o ${my_object_files} \ + ${LAMMPS_LIB} ${MPI_STUBS} -lstdc++ -lm + OR + ifort LAMMPS.o LAMMPS-wrapper.o ${my_object_files} \ + ${LAMMPS_LIB} ${MPI_STUBS} -cxxlib -lm + Any other required libraries (e.g. -ljpeg, -lfftw) should be added to + the end of this line. + +You should now have a working executable. + +** DYNAMIC LIBRARY INSTRUCTIONS ** + (1) Compile LAMMPS as a dynamic library + (make makeshlib && make -f Makefile.shlib [targetname]). + (2) Compile, but don't link, LAMMPS.F90 using the -fPIC flag, such as + mpif90 -fPIC -c LAMMPS.f90 + (3) Compile, but don't link, LAMMPS-wrapper.cpp in the same manner, e.g. + mpicxx -fPIC -c LAMMPS-wrapper.cpp + (4) Make the dynamic library, like so: + mpif90 -fPIC -shared -o liblammps_fortran.so LAMMPS.o LAMMPS-wrapper.o + (5) Compile your program, such as, + mpif90 -I${LAMMPS_SRC}/examples/COUPLE/fortran2 -c myfreeformatfile.f90 + where ${LAMMPS_SRC}/examples/COUPLE/fortran2 contains the .mod file from + step (3) + (6) Link everything together, such as + mpif90 ${my_object_files} -L${LAMMPS_SRC} \ + -L${LAMMPS_SRC}/examples/COUPLE/fortran2 -llammps_fortran \ + -llammps_openmpi -lmpi_cxx -lstdc++ -lm + +If you wish to avoid the -L flags, add the directories containing your +shared libraries to the LIBRARY_PATH environment variable. At run time, you +will have to add these directories to LD_LIBRARY_PATH as well; otherwise, +your executable will not find the libraries it needs. + +------------------------------------- + +--USAGE-- + +To use this API, your program unit (PROGRAM/SUBROUTINE/FUNCTION/MODULE/etc.) +should look something like this: + program call_lammps + use LAMMPS + ! Other modules, etc. + implicit none + type (lammps_instance) :: lmp ! This is a pointer to your LAMMPS instance + real (C_double) :: fix + real (C_double), dimension(:), pointer :: fix2 + ! Rest of declarations + call lammps_open_no_mpi ('lmp -in /dev/null -screen out.lammps',lmp) + ! Set up rest of program here + call lammps_file (lmp, 'in.example') + call lammps_extract_fix (fix, lmp, '2', 0, 1, 1, 1) + call lammps_extract_fix (fix2, lmp, '4', 0, 2, 1, 1) + call lammps_close (lmp) + end program call_lammps + +Important notes: + * Though I dislike the use of pointers, they are necessary when communicating + with C and C++, which do not support Fortran's ALLOCATABLE attribute. + * There is no need to deallocate C-allocated memory; this is done for you in + the cases when it is done (which are all cases when pointers are not + accepted, such as global fix data) + * All arguments which are char* variables in library.cpp are character (len=*) + variables here. For example, + call lammps_command (lmp, 'units metal') + will work as expected. + * The public functions (the only ones you can use) have interfaces as + described in the comments at the top of LAMMPS.F90. They are not always + the same as those in library.h, since C strings are replaced by Fortran + strings and the like. + * The module attempts to check whether you have done something stupid (such + as assign a 2D array to a scalar), but it's not perfect. For example, the + command + call lammps_extract_global (nlocal, ptr, 'nlocal') + will give nlocal correctly if nlocal is a pointer to type INTEGER, but it + will give the wrong answer if nlocal is a pointer to type REAL. This is a + feature of the (void*) type cast in library.cpp. There is no way I can + check this for you! It WILL catch you if you pass it an allocatable or + fixed-size array when it expects a pointer. + * Arrays constructed from temporary data from LAMMPS are ALLOCATABLE, and + represent COPIES of data, not the originals. Functions like + lammps_extract_atom, which return actual LAMMPS data, are pointers. + * IMPORTANT: Due to the differences between C and Fortran arrays (C uses + row-major vectors, Fortran uses column-major vectors), all arrays returned + from LAMMPS have their indices swapped. + * An example of a complete program, simple.f90, is included with this + package. + +------------------------------------- + +--TROUBLESHOOTING-- + +Compile-time errors (when compiling LAMMPS.F90, that is) probably indicate +that your compiler is not new enough to support Fortran 2003 features. For +example, GCC 4.1.2 will not compile this module, but GCC 4.4.0 will. + +If your compiler balks at 'use, intrinsic :: ISO_C_binding,' try removing the +intrinsic part so it looks like an ordinary module. However, it is likely +that such a compiler will also have problems with everything else in the +file as well. + +If you get a segfault as soon as the lammps_open call is made, check that you +compiled your program AND LAMMPS-wrapper.cpp using the same MPI headers. Using +the stubs for one and the actual MPI library for the other will cause Bad +Things to happen. + +If you find run-time errors, please pass them along via the LAMMPS Users +mailing list (please CC me as well; address above). Please provide a minimal +working example along with the names and versions of the compilers you are +using. Please make sure the error is repeatable and is in MY code, not yours +(generating a minimal working example will usually ensure this anyway). + +------------------------------------- + +--IMPLEMENTATION NOTES-- + +The Fortran procedures have the same names as the C procedures, and +their purpose is the same, but they may take different arguments. Here are +some of the important differences: + * lammps_open and lammps_open_no_mpi take a string instead of argc and + argv. This is necessary because C and C++ have a very different way + of treating strings than Fortran. If you want the command line to be + passed to lammps_open (as it often would be from C/C++), use the + GET_COMMAND intrinsic to obtain it. + * All C++ functions that accept char* pointers now accept Fortran-style + strings within this interface instead. + * All of the lammps_extract_[something] functions, which return void* + C-style pointers, have been replaced by generic subroutines that return + Fortran variables (which may be arrays). The first argument houses the + variable/pointer to be returned (pretend it's on the left-hand side); all + other arguments are identical except as stipulated above. + Note that it is not possible to declare generic functions that are selected + based solely on the type/kind/rank (TKR) signature of the return value, + only based on the TKR of the arguments. + * The SHAPE of the first argument to lammps_extract_[something] is checked + against the "shape" of the C array (e.g., double vs. double* vs. double**). + Calling a subroutine with arguments of inappropriate rank will result in an + error at run time. + * The indices i and j in lammps_extract_fix are used the same way they + are in f_ID[i][j] references in LAMMPS (i.e., starting from 1). This is + different than the way library.cpp uses these numbers, but is more + consistent with the way arrays are accessed in LAMMPS and in Fortran. + * The char* pointer normally returned by lammps_command is thrown away + in this version; note also that lammps_command is now a subroutine + instead of a function. + * The pointer to LAMMPS itself is of type(lammps_instance), which is itself + a synonym for type(C_ptr), part of ISO_C_BINDING. Type (C_ptr) is + C's void* data type. + * This module will almost certainly generate a compile-time warning, + such as, + subroutine lammps_open_wrapper (argc, argv, communicator, ptr) & + Variable 'communicator' at (1) is a parameter to the BIND(C) + procedure 'lammps_open_wrapper' but may not be C interoperable + This happens because lammps_open_wrapper actually takes a Fortran + INTEGER argument, whose type is defined by the MPI library itself. The + Fortran integer is converted to a C integer by the MPI library (if such + conversion is actually necessary). + * lammps_extract_global returns COPIES of the (scalar) data, as does the + C version. + * lammps_extract_atom, lammps_extract_compute, and lammps_extract_fix + have a first argument that will be associated with ACTUAL LAMMPS DATA. + This means the first argument must be: + * The right rank (via the DIMENSION modifier) + * A C-interoperable POINTER type (i.e., INTEGER (C_int) or + REAL (C_double)). + * lammps_extract_variable returns COPIES of the data, as the C library + interface does. There is no need to deallocate using lammps_free. + * The 'data' argument to lammps_gather_atoms and lammps_scatter atoms must + be ALLOCATABLE. It should be of type INTEGER or DOUBLE PRECISION. It + does NOT need to be C inter-operable (and indeed should not be). + * The 'count' argument of lammps_scatter_atoms is unnecessary; the shape of + the array determines the number of elements LAMMPS will read. diff --git a/examples/COUPLE/fortran2/in.simple b/examples/COUPLE/fortran2/in.simple index 69384836ee..5982cbaac1 100644 --- a/examples/COUPLE/fortran2/in.simple +++ b/examples/COUPLE/fortran2/in.simple @@ -1,10 +1,11 @@ -units metal -lattice bcc 3.1656 +units lj +atom_modify map array +lattice bcc 1.0 region simbox block 0 10 0 10 0 10 create_box 2 simbox create_atoms 1 region simbox -pair_style eam/fs -pair_coeff * * path/to/my_potential.eam.fs A1 A2 +pair_style lj/cut 2.5 +pair_coeff * * 1.0 1.0 mass 1 58.2 # These are made-up numbers mass 2 28.3 velocity all create 1200.0 7474848 dist gaussian diff --git a/examples/COUPLE/fortran2/simple.f90 b/examples/COUPLE/fortran2/simple.f90 index 93cf519f97..7ed3850a3d 100644 --- a/examples/COUPLE/fortran2/simple.f90 +++ b/examples/COUPLE/fortran2/simple.f90 @@ -1,44 +1,111 @@ program simple + use MPI use LAMMPS + + ! The following line is unnecessary, as I have included these three entities + ! with the LAMMPS module, but I leave them in anyway to remind people where + ! they came from + use, intrinsic :: ISO_C_binding, only : C_double, C_ptr, C_int + implicit none - type (lammps_instance) :: lmp - double precision :: compute, fix, fix2 - double precision, dimension(:), allocatable :: compute_v, mass, r - double precision, dimension(:,:), allocatable :: x - real, dimension(:,:), allocatable :: x_r + ! Notes: + ! * If LAMMPS returns a scalar that is allocated by the library interface + ! (see library.cpp), then that memory is deallocated automatically and + ! the argument to lammps_extract_fix must be a SCALAR. + ! * If LAMMPS returns a pointer to an array, consisting of internal LAMMPS + ! data, then the argument must be an interoperable Fortran pointer. + ! Interoperable means it is of type INTEGER (C_INT) or of type + ! REAL (C_DOUBLE) in this context. + ! * Pointers should NEVER be deallocated, as that would deallocate internal + ! LAMMPS data! + ! * Note that just because you can read the values of, say, a compute at + ! any time does not mean those values represent the "correct" values. + ! LAMMPS will abort you if you try to grab a pointer to a non-current + ! entity, but once it's bound, it's your responsibility to check that + ! it's current before evaluating. + ! * IMPORTANT: Two-dimensional arrays (such as 'x' from extract_atom) + ! will be transposed from what they might look like in C++. This is + ! because of different bookkeeping conventions between Fortran and C + ! that date back to about 1970 or so (when C was written). + ! * Arrays start from 1, EXCEPT for mass from extract_atom, which + ! starts from 0. This is because the C array actually has a blank + ! first element (and thus mass[1] corresponds to the mass of type 1) + + type (C_ptr) :: lmp + real (C_double), pointer :: compute => NULL() + real (C_double) :: fix, fix2 + real (C_double), dimension(:), pointer :: compute_v => NULL() + real (C_double), dimension(:,:), pointer :: x => NULL() + real (C_double), dimension(:), pointer :: mass => NULL() + integer, dimension(:), allocatable :: types + double precision, dimension(:), allocatable :: r + integer :: error, narg, me, nprocs + character (len=1024) :: command_line + + call MPI_Init (error) + call MPI_Comm_rank (MPI_COMM_WORLD, me, error) + call MPI_Comm_size (MPI_COMM_WORLD, nprocs, error) + + ! You are free to pass any string you like to lammps_open or + ! lammps_open_no_mpi; here is how you pass it the command line + !call get_command (command_line) + !call lammps_open (command_line, MPI_COMM_WORLD, lmp) + + ! And here's how to to it with a string constant of your choice + call lammps_open_no_mpi ('lmp -log log.simple', lmp) - call lammps_open_no_mpi ('',lmp) call lammps_file (lmp, 'in.simple') call lammps_command (lmp, 'run 500') + ! This extracts f_2 as a scalar (the last two arguments can be arbitrary) call lammps_extract_fix (fix, lmp, '2', 0, 1, 1, 1) print *, 'Fix is ', fix + ! This extracts f_4[1][1] as a scalar call lammps_extract_fix (fix2, lmp, '4', 0, 2, 1, 1) print *, 'Fix 2 is ', fix2 + ! This extracts the scalar compute of compute thermo_temp call lammps_extract_compute (compute, lmp, 'thermo_temp', 0, 0) print *, 'Compute is ', compute + ! This extracts the vector compute of compute thermo_temp call lammps_extract_compute (compute_v, lmp, 'thermo_temp', 0, 1) print *, 'Vector is ', compute_v + ! This extracts the masses call lammps_extract_atom (mass, lmp, 'mass') - print *, 'Mass is ', mass + print *, 'Mass is ', mass(1:) + ! Extracts a pointer to the arrays of positions for all atoms call lammps_extract_atom (x, lmp, 'x') - if ( .not. allocated (x) ) print *, 'x is not allocated' - print *, 'x is ', x(1,:) + if ( .not. associated (x) ) print *, 'x is not associated' + print *, 'x is ', x(:,1) ! Prints x, y, z for atom 1 - call lammps_extract_atom (x_r, lmp, 'x') - if ( .not. allocated (x_r) ) print *, 'x is not allocated' - print *, 'x_r is ', x_r(1,:) + ! Extracts pointer to atom types + call lammps_gather_atoms (lmp, 'type', 1, types) + print *, 'types is ', types(1:3) - call lammps_get_coords (lmp, r) - print *, 'r is ', r(1:3) + ! Allocates an array and assigns all positions to it + call lammps_gather_atoms (lmp, 'x', 3, r) + print *, 'size(r) = ', size(r) + print *, 'r is ', r(1:6) + + ! Puts those position data back + call lammps_scatter_atoms (lmp, 'x', r) + + call lammps_command (lmp, 'run 1') + print *, 'x is ', x(:,1) ! Note that the position updates! + print *, 'Compute is ', compute ! This did only because "temp" is part of + ! the thermo output; the vector part did + ! not, and won't until we give LAMMPS a + ! thermo output or other command that + ! requires its value call lammps_close (lmp) + call MPI_Finalize (error) + end program simple From cb9d45386c198e5de3ee2f4dad666b0de8cdf454 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 25 Oct 2012 15:47:22 +0000 Subject: [PATCH 06/12] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9004 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/version.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/version.h b/src/version.h index 7cfd1a5c93..5d2c676261 100644 --- a/src/version.h +++ b/src/version.h @@ -1 +1 @@ -#define LAMMPS_VERSION "23 Oct 2012" +#define LAMMPS_VERSION "25 Oct 2012" From 6f4bba7edf244656cd1a552d8765c499646acc73 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 25 Oct 2012 16:24:37 +0000 Subject: [PATCH 07/12] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9006 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/compute_com_molecule.cpp | 16 +-- src/compute_gyration.cpp | 17 ++- src/compute_gyration_molecule.cpp | 57 ++++------ src/compute_msd_molecule.cpp | 22 ++-- src/fix_addforce.cpp | 15 +-- src/fix_momentum.cpp | 18 ++- src/fix_spring_rg.cpp | 15 +-- src/fix_spring_self.cpp | 34 ++---- src/fix_tmd.cpp | 51 +++------ src/group.cpp | 176 +++++++++++------------------- src/velocity.cpp | 15 +-- 11 files changed, 149 insertions(+), 287 deletions(-) diff --git a/src/compute_com_molecule.cpp b/src/compute_com_molecule.cpp index 80e745f056..82a212ca3a 100644 --- a/src/compute_com_molecule.cpp +++ b/src/compute_com_molecule.cpp @@ -96,8 +96,8 @@ void ComputeCOMMolecule::init() void ComputeCOMMolecule::compute_array() { int i,imol; - double xbox,ybox,zbox; double massone; + double unwrap[3]; invoked_array = update->ntimestep; @@ -113,23 +113,17 @@ void ComputeCOMMolecule::compute_array() double *rmass = atom->rmass; int nlocal = atom->nlocal; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; imol = molecule[i]; if (molmap) imol = molmap[imol-idlo]; else imol--; - com[imol][0] += (x[i][0] + xbox*xprd) * massone; - com[imol][1] += (x[i][1] + ybox*yprd) * massone; - com[imol][2] += (x[i][2] + zbox*zprd) * massone; + domain->unmap(x[i],image[i],unwrap); + com[imol][0] += unwrap[0] * massone; + com[imol][1] += unwrap[1] * massone; + com[imol][2] += unwrap[2] * massone; } MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules, diff --git a/src/compute_gyration.cpp b/src/compute_gyration.cpp index 2e65cf1793..007d56f0c3 100644 --- a/src/compute_gyration.cpp +++ b/src/compute_gyration.cpp @@ -83,25 +83,22 @@ void ComputeGyration::compute_vector() double *rmass = atom->rmass; int nlocal = atom->nlocal; - int xbox,ybox,zbox; double dx,dy,dz,massone; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; + double unwrap[3]; double rg[6]; rg[0] = rg[1] = rg[2] = rg[3] = rg[4] = rg[5] = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - xcm[0]; - dy = (x[i][1] + ybox*yprd) - xcm[1]; - dz = (x[i][2] + zbox*zprd) - xcm[2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; + + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - xcm[0]; + dy = unwrap[1] - xcm[1]; + dz = unwrap[2] - xcm[2]; + rg[0] += dx*dx * massone; rg[1] += dy*dy * massone; rg[2] += dz*dz * massone; diff --git a/src/compute_gyration_molecule.cpp b/src/compute_gyration_molecule.cpp index 8fd9b6751e..1ebc47a9aa 100644 --- a/src/compute_gyration_molecule.cpp +++ b/src/compute_gyration_molecule.cpp @@ -123,8 +123,8 @@ void ComputeGyrationMolecule::init() void ComputeGyrationMolecule::compute_vector() { int i,imol; - double xbox,ybox,zbox,dx,dy,dz; - double massone; + double dx,dy,dz,massone; + double unwrap[3]; invoked_array = update->ntimestep; @@ -141,21 +141,16 @@ void ComputeGyrationMolecule::compute_vector() double *rmass = atom->rmass; int nlocal = atom->nlocal; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; imol = molecule[i]; if (molmap) imol = molmap[imol-idlo]; else imol--; - dx = (x[i][0] + xbox*xprd) - comall[imol][0]; - dy = (x[i][1] + ybox*yprd) - comall[imol][1]; - dz = (x[i][2] + zbox*zprd) - comall[imol][2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - comall[imol][0]; + dy = unwrap[1] - comall[imol][1]; + dz = unwrap[2] - comall[imol][2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; rg[imol] += (dx*dx + dy*dy + dz*dz) * massone; @@ -171,8 +166,8 @@ void ComputeGyrationMolecule::compute_vector() void ComputeGyrationMolecule::compute_array() { int i,j,imol; - double xbox,ybox,zbox,dx,dy,dz; - double massone; + double dx,dy,dz,massone; + double unwrap[3]; invoked_array = update->ntimestep; @@ -191,21 +186,15 @@ void ComputeGyrationMolecule::compute_array() double *rmass = atom->rmass; int nlocal = atom->nlocal; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; imol = molecule[i]; if (molmap) imol = molmap[imol-idlo]; else imol--; - dx = (x[i][0] + xbox*xprd) - comall[imol][0]; - dy = (x[i][1] + ybox*yprd) - comall[imol][1]; - dz = (x[i][2] + zbox*zprd) - comall[imol][2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - comall[imol][0]; + dy = unwrap[1] - comall[imol][1]; + dz = unwrap[2] - comall[imol][2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; rgt[imol][0] += dx*dx * massone; @@ -233,8 +222,8 @@ void ComputeGyrationMolecule::compute_array() void ComputeGyrationMolecule::molcom() { int i,imol; - double xbox,ybox,zbox,dx,dy,dz; - double massone; + double dx,dy,dz,massone; + double unwrap[3]; for (i = 0; i < nmolecules; i++) com[i][0] = com[i][1] = com[i][2] = 0.0; @@ -248,23 +237,17 @@ void ComputeGyrationMolecule::molcom() double *rmass = atom->rmass; int nlocal = atom->nlocal; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; imol = molecule[i]; if (molmap) imol = molmap[imol-idlo]; else imol--; - com[imol][0] += (x[i][0] + xbox*xprd) * massone; - com[imol][1] += (x[i][1] + ybox*yprd) * massone; - com[imol][2] += (x[i][2] + zbox*zprd) * massone; + domain->unmap(x[i],image[i],unwrap); + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + com[imol][0] += unwrap[0] * massone; + com[imol][1] += unwrap[1] * massone; + com[imol][2] += unwrap[2] * massone; } MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules, diff --git a/src/compute_msd_molecule.cpp b/src/compute_msd_molecule.cpp index c98cca8e1e..270cfad6f8 100644 --- a/src/compute_msd_molecule.cpp +++ b/src/compute_msd_molecule.cpp @@ -111,8 +111,8 @@ void ComputeMSDMolecule::init() void ComputeMSDMolecule::compute_array() { int i,imol; - double xbox,ybox,zbox,dx,dy,dz; - double massone; + double dx,dy,dz,massone; + double unwrap[3]; invoked_array = update->ntimestep; @@ -130,23 +130,17 @@ void ComputeMSDMolecule::compute_array() double *rmass = atom->rmass; int nlocal = atom->nlocal; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; imol = molecule[i]; if (molmap) imol = molmap[imol-idlo]; else imol--; - com[imol][0] += (x[i][0] + xbox*xprd) * massone; - com[imol][1] += (x[i][1] + ybox*yprd) * massone; - com[imol][2] += (x[i][2] + zbox*zprd) * massone; + domain->unmap(x[i],image[i],unwrap); + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + com[imol][0] += unwrap[0] * massone; + com[imol][1] += unwrap[1] * massone; + com[imol][2] += unwrap[2] * massone; } MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules, diff --git a/src/fix_addforce.cpp b/src/fix_addforce.cpp index 2523c7cc05..8ec546bb01 100644 --- a/src/fix_addforce.cpp +++ b/src/fix_addforce.cpp @@ -217,12 +217,6 @@ void FixAddForce::min_setup(int vflag) void FixAddForce::post_force(int vflag) { - int xbox,ybox,zbox; - - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - double **x = atom->x; double **f = atom->f; int *mask = atom->mask; @@ -247,18 +241,15 @@ void FixAddForce::post_force(int vflag) // potential energy = - x dot f in unwrapped coords if (varflag == CONSTANT) { + double unwrap[3]; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (iregion >= 0 && !domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) continue; - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - foriginal[0] -= xvalue * (x[i][0]+xbox*xprd) + - yvalue * (x[i][1]+ybox*yprd) + zvalue * (x[i][2]+zbox*zprd); - + domain->unmap(x[i],image[i],unwrap); + foriginal[0] -= xvalue*unwrap[0] + yvalue*unwrap[1] + zvalue*unwrap[2]; foriginal[1] += f[i][0]; foriginal[2] += f[i][1]; foriginal[3] += f[i][2]; diff --git a/src/fix_momentum.cpp b/src/fix_momentum.cpp index 69e06aeaed..6f9da4427b 100644 --- a/src/fix_momentum.cpp +++ b/src/fix_momentum.cpp @@ -57,7 +57,8 @@ FixMomentum::FixMomentum(LAMMPS *lmp, int narg, char **arg) : if (linear) if (xflag < 0 || xflag > 1 || yflag < 0 || yflag > 1 || - zflag < 0 || zflag > 1) error->all(FLERR,"Illegal fix momentum command"); + zflag < 0 || zflag > 1) + error->all(FLERR,"Illegal fix momentum command"); // cannot have 0 atoms in group @@ -121,20 +122,15 @@ void FixMomentum::end_of_step() tagint *image = atom->image; int nlocal = atom->nlocal; - int xbox,ybox,zbox; double dx,dy,dz; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; + double unwrap[3]; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - xcm[0]; - dy = (x[i][1] + ybox*yprd) - xcm[1]; - dz = (x[i][2] + zbox*zprd) - xcm[2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - xcm[0]; + dy = unwrap[1] - xcm[1]; + dz = unwrap[2] - xcm[2]; v[i][0] -= omega[1]*dz - omega[2]*dy; v[i][1] -= omega[2]*dx - omega[0]*dz; v[i][2] -= omega[0]*dy - omega[1]*dx; diff --git a/src/fix_spring_rg.cpp b/src/fix_spring_rg.cpp index e17aa8966c..0d93c4b628 100644 --- a/src/fix_spring_rg.cpp +++ b/src/fix_spring_rg.cpp @@ -110,20 +110,15 @@ void FixSpringRG::post_force(int vflag) int nlocal = atom->nlocal; double massfrac; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; + double unwrap[3]; - int xbox,ybox,zbox; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - xcm[0]; + dy = unwrap[1] - xcm[1]; + dz = unwrap[2] - xcm[2]; term1 = 2.0 * k * (1.0 - rg0/rg); - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - xcm[0]; - dy = (x[i][1] + ybox*yprd) - xcm[1]; - dz = (x[i][2] + zbox*zprd) - xcm[2]; massfrac = mass[type[i]]/masstotal; f[i][0] -= term1*dx*massfrac; f[i][1] -= term1*dy*massfrac; diff --git a/src/fix_spring_self.cpp b/src/fix_spring_self.cpp index fb2120d4cb..043bd21e2e 100644 --- a/src/fix_spring_self.cpp +++ b/src/fix_spring_self.cpp @@ -45,9 +45,10 @@ FixSpringSelf::FixSpringSelf(LAMMPS *lmp, int narg, char **arg) : if (k <= 0.0) error->all(FLERR,"Illegal fix spring/self command"); xflag = yflag = zflag = 1; + if (narg == 5) { if (strcmp(arg[4],"xyz") == 0) { - ; /* default */ + xflag = yflag = zflag = 1; } else if (strcmp(arg[4],"xy") == 0) { zflag = 0; } else if (strcmp(arg[4],"xz") == 0) { @@ -78,20 +79,9 @@ FixSpringSelf::FixSpringSelf(LAMMPS *lmp, int narg, char **arg) : tagint *image = atom->image; int nlocal = atom->nlocal; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - int xbox,ybox,zbox; - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - xoriginal[i][0] = x[i][0] + xbox*xprd; - xoriginal[i][1] = x[i][1] + ybox*yprd; - xoriginal[i][2] = x[i][2] + zbox*zprd; - } else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0; + if (mask[i] & groupbit) domain->unmap(x[i],image[i],xoriginal[i]); + else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0; } espring = 0.0; @@ -161,21 +151,17 @@ void FixSpringSelf::post_force(int vflag) tagint *image = atom->image; int nlocal = atom->nlocal; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - int xbox,ybox,zbox; double dx,dy,dz; + double unwrap[3]; + espring = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = x[i][0] + xbox*xprd - xoriginal[i][0]; - dy = x[i][1] + ybox*yprd - xoriginal[i][1]; - dz = x[i][2] + zbox*zprd - xoriginal[i][2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - xoriginal[i][0]; + dy = unwrap[1] - xoriginal[i][1]; + dz = unwrap[2] - xoriginal[i][2]; if (!xflag) dx = 0.0; if (!yflag) dy = 0.0; if (!zflag) dz = 0.0; diff --git a/src/fix_tmd.cpp b/src/fix_tmd.cpp index 71efd01bdf..fde131e632 100644 --- a/src/fix_tmd.cpp +++ b/src/fix_tmd.cpp @@ -98,26 +98,17 @@ FixTMD::FixTMD(LAMMPS *lmp, int narg, char **arg) : double *mass = atom->mass; int nlocal = atom->nlocal; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - double dx,dy,dz; - int xbox,ybox,zbox; + rho_start = 0.0; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = x[i][0] + xbox*xprd - xf[i][0]; - dy = x[i][1] + ybox*yprd - xf[i][1]; - dz = x[i][2] + zbox*zprd - xf[i][2]; + domain->unmap(x[i],image[i],xold[i]); + dx = xold[i][0] - xf[i][0]; + dy = xold[i][1] - xf[i][1]; + dz = xold[i][2] - xf[i][2]; rho_start += mass[type[i]]*(dx*dx + dy*dy + dz*dz); - xold[i][0] = x[i][0] + xbox*xprd; - xold[i][1] = x[i][1] + ybox*yprd; - xold[i][2] = x[i][2] + zbox*zprd; } else xold[i][0] = xold[i][1] = xold[i][2] = 0.0; } @@ -190,10 +181,8 @@ void FixTMD::initial_integrate(int vflag) double dxold,dyold,dzold,xback,yback,zback; double gamma_forward,gamma_back,gamma_max,lambda; double kt,fr,kttotal,frtotal,dtfm; + double unwrap[3]; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; double **x = atom->x; double **v = atom->v; double **f = atom->f; @@ -202,7 +191,6 @@ void FixTMD::initial_integrate(int vflag) int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; - int xbox,ybox,zbox; double delta = update->ntimestep - update->beginstep; delta /= update->endstep - update->beginstep; @@ -216,12 +204,10 @@ void FixTMD::initial_integrate(int vflag) dxold = xold[i][0] - xf[i][0]; dyold = xold[i][1] - xf[i][1]; dzold = xold[i][2] - xf[i][2]; - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = x[i][0] + xbox*xprd - xf[i][0]; - dy = x[i][1] + ybox*yprd - xf[i][1]; - dz = x[i][2] + zbox*zprd - xf[i][2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - xf[i][0]; + dy = unwrap[1] - xf[i][1]; + dz = unwrap[2] - xf[i][2]; a += mass[type[i]]*(dxold*dxold + dyold*dyold + dzold*dzold); b += mass[type[i]]*(dx *dxold + dy *dyold + dz *dzold); e += mass[type[i]]*(dx *dx + dy *dy + dz *dz); @@ -258,12 +244,10 @@ void FixTMD::initial_integrate(int vflag) dxold = xold[i][0] - xf[i][0]; dyold = xold[i][1] - xf[i][1]; dzold = xold[i][2] - xf[i][2]; - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - xback = x[i][0] + xbox*xprd + gamma_back*dxold; - yback = x[i][1] + ybox*yprd + gamma_back*dyold; - zback = x[i][2] + zbox*zprd + gamma_back*dzold; + domain->unmap(x[i],image[i],unwrap); + xback = unwrap[0] + gamma_back*dxold; + yback = unwrap[1] + gamma_back*dyold; + zback = unwrap[2] + gamma_back*dzold; dxkt = xback - xold[i][0]; dykt = yback - xold[i][1]; dzkt = zback - xold[i][2]; @@ -314,12 +298,7 @@ void FixTMD::initial_integrate(int vflag) x[i][2] += gamma_forward*dzold; v[i][2] += gamma_forward*dzold/dtv; f[i][2] += gamma_forward*dzold/dtv/dtfm; - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - xold[i][0] = x[i][0] + xbox*xprd; - xold[i][1] = x[i][1] + ybox*yprd; - xold[i][2] = x[i][2] + zbox*zprd; + domain->unmap(x[i],image[i],xold[i]); } } } diff --git a/src/group.cpp b/src/group.cpp index 21bea8f750..c82d110058 100644 --- a/src/group.cpp +++ b/src/group.cpp @@ -766,34 +766,27 @@ void Group::xcm(int igroup, double masstotal, double *cm) double cmone[3]; cmone[0] = cmone[1] = cmone[2] = 0.0; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - int xbox,ybox,zbox; double massone; + double unwrap[3]; if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; massone = rmass[i]; - cmone[0] += (x[i][0] + xbox*xprd) * massone; - cmone[1] += (x[i][1] + ybox*yprd) * massone; - cmone[2] += (x[i][2] + zbox*zprd) * massone; + domain->unmap(x[i],image[i],unwrap); + cmone[0] += unwrap[0] * massone; + cmone[1] += unwrap[1] * massone; + cmone[2] += unwrap[2] * massone; } } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; massone = mass[type[i]]; - cmone[0] += (x[i][0] + xbox*xprd) * massone; - cmone[1] += (x[i][1] + ybox*yprd) * massone; - cmone[2] += (x[i][2] + zbox*zprd) * massone; + domain->unmap(x[i],image[i],unwrap); + cmone[0] += unwrap[0] * massone; + cmone[1] += unwrap[1] * massone; + cmone[2] += unwrap[2] * massone; } } @@ -827,34 +820,27 @@ void Group::xcm(int igroup, double masstotal, double *cm, int iregion) double cmone[3]; cmone[0] = cmone[1] = cmone[2] = 0.0; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - int xbox,ybox,zbox; double massone; + double unwrap[3]; if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; massone = rmass[i]; - cmone[0] += (x[i][0] + xbox*xprd) * massone; - cmone[1] += (x[i][1] + ybox*yprd) * massone; - cmone[2] += (x[i][2] + zbox*zprd) * massone; + domain->unmap(x[i],image[i],unwrap); + cmone[0] += unwrap[0] * massone; + cmone[1] += unwrap[1] * massone; + cmone[2] += unwrap[2] * massone; } } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; massone = mass[type[i]]; - cmone[0] += (x[i][0] + xbox*xprd) * massone; - cmone[1] += (x[i][1] + ybox*yprd) * massone; - cmone[2] += (x[i][2] + zbox*zprd) * massone; + domain->unmap(x[i],image[i],unwrap); + cmone[0] += unwrap[0] * massone; + cmone[1] += unwrap[1] * massone; + cmone[2] += unwrap[2] * massone; } } @@ -1102,21 +1088,16 @@ double Group::gyration(int igroup, double masstotal, double *cm) double *rmass = atom->rmass; int nlocal = atom->nlocal; - int xbox,ybox,zbox; double dx,dy,dz,massone; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; + double unwrap[3]; double rg = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - cm[0]; - dy = (x[i][1] + ybox*yprd) - cm[1]; - dz = (x[i][2] + zbox*zprd) - cm[2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - cm[0]; + dy = unwrap[1] - cm[1]; + dz = unwrap[2] - cm[2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; rg += (dx*dx + dy*dy + dz*dz) * massone; @@ -1147,21 +1128,16 @@ double Group::gyration(int igroup, double masstotal, double *cm, int iregion) double *rmass = atom->rmass; int nlocal = atom->nlocal; - int xbox,ybox,zbox; double dx,dy,dz,massone; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; + double unwrap[3]; double rg = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - cm[0]; - dy = (x[i][1] + ybox*yprd) - cm[1]; - dz = (x[i][2] + zbox*zprd) - cm[2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - cm[0]; + dy = unwrap[1] - cm[1]; + dz = unwrap[2] - cm[2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; rg += (dx*dx + dy*dy + dz*dz) * massone; @@ -1192,22 +1168,18 @@ void Group::angmom(int igroup, double *cm, double *lmom) double *rmass = atom->rmass; int nlocal = atom->nlocal; - int xbox,ybox,zbox; double dx,dy,dz,massone; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; + double unwrap[3]; + double p[3]; p[0] = p[1] = p[2] = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - cm[0]; - dy = (x[i][1] + ybox*yprd) - cm[1]; - dz = (x[i][2] + zbox*zprd) - cm[2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - cm[0]; + dy = unwrap[1] - cm[1]; + dz = unwrap[2] - cm[2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; p[0] += massone * (dy*v[i][2] - dz*v[i][1]); @@ -1238,22 +1210,18 @@ void Group::angmom(int igroup, double *cm, double *lmom, int iregion) double *rmass = atom->rmass; int nlocal = atom->nlocal; - int xbox,ybox,zbox; double dx,dy,dz,massone; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; + double unwrap[3]; + double p[3]; p[0] = p[1] = p[2] = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - cm[0]; - dy = (x[i][1] + ybox*yprd) - cm[1]; - dz = (x[i][2] + zbox*zprd) - cm[2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - cm[0]; + dy = unwrap[1] - cm[1]; + dz = unwrap[2] - cm[2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; p[0] += massone * (dy*v[i][2] - dz*v[i][1]); @@ -1280,22 +1248,18 @@ void Group::torque(int igroup, double *cm, double *tq) tagint *image = atom->image; int nlocal = atom->nlocal; - int xbox,ybox,zbox; double dx,dy,dz; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; + double unwrap[3]; + double tlocal[3]; tlocal[0] = tlocal[1] = tlocal[2] = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - cm[0]; - dy = (x[i][1] + ybox*yprd) - cm[1]; - dz = (x[i][2] + zbox*zprd) - cm[2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - cm[0]; + dy = unwrap[1] - cm[1]; + dz = unwrap[2] - cm[2]; tlocal[0] += dy*f[i][2] - dz*f[i][1]; tlocal[1] += dz*f[i][0] - dx*f[i][2]; tlocal[2] += dx*f[i][1] - dy*f[i][0]; @@ -1321,22 +1285,18 @@ void Group::torque(int igroup, double *cm, double *tq, int iregion) tagint *image = atom->image; int nlocal = atom->nlocal; - int xbox,ybox,zbox; double dx,dy,dz; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; + double unwrap[3]; + double tlocal[3]; tlocal[0] = tlocal[1] = tlocal[2] = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - cm[0]; - dy = (x[i][1] + ybox*yprd) - cm[1]; - dz = (x[i][2] + zbox*zprd) - cm[2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - cm[0]; + dy = unwrap[1] - cm[1]; + dz = unwrap[2] - cm[2]; tlocal[0] += dy*f[i][2] - dz*f[i][1]; tlocal[1] += dz*f[i][0] - dx*f[i][2]; tlocal[2] += dx*f[i][1] - dy*f[i][0]; @@ -1364,11 +1324,9 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3]) double *rmass = atom->rmass; int nlocal = atom->nlocal; - int xbox,ybox,zbox; double dx,dy,dz,massone; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; + double unwrap[3]; + double ione[3][3]; for (i = 0; i < 3; i++) for (j = 0; j < 3; j++) @@ -1376,12 +1334,10 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3]) for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - cm[0]; - dy = (x[i][1] + ybox*yprd) - cm[1]; - dz = (x[i][2] + zbox*zprd) - cm[2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - cm[0]; + dy = unwrap[1] - cm[1]; + dz = unwrap[2] - cm[2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; ione[0][0] += massone * (dy*dy + dz*dz); @@ -1418,11 +1374,9 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3], int iregion) double *rmass = atom->rmass; int nlocal = atom->nlocal; - int xbox,ybox,zbox; double dx,dy,dz,massone; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; + double unwrap[3]; + double ione[3][3]; for (i = 0; i < 3; i++) for (j = 0; j < 3; j++) @@ -1430,12 +1384,10 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3], int iregion) for (i = 0; i < nlocal; i++) if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - cm[0]; - dy = (x[i][1] + ybox*yprd) - cm[1]; - dz = (x[i][2] + zbox*zprd) - cm[2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - cm[0]; + dy = unwrap[1] - cm[1]; + dz = unwrap[2] - cm[2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; ione[0][0] += massone * (dy*dy + dz*dz); diff --git a/src/velocity.cpp b/src/velocity.cpp index 6063eb46fc..e69b6d67a5 100644 --- a/src/velocity.cpp +++ b/src/velocity.cpp @@ -718,20 +718,15 @@ void Velocity::zero_rotation() tagint *image = atom->image; int nlocal = atom->nlocal; - int xbox,ybox,zbox; double dx,dy,dz; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; + double unwrap[3]; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - xcm[0]; - dy = (x[i][1] + ybox*yprd) - xcm[1]; - dz = (x[i][2] + zbox*zprd) - xcm[2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0]- xcm[0]; + dy = unwrap[1] - xcm[1]; + dz = unwrap[2] - xcm[2]; v[i][0] -= omega[1]*dz - omega[2]*dy; v[i][1] -= omega[2]*dx - omega[0]*dz; v[i][2] -= omega[0]*dy - omega[1]*dx; From 43b161eb3398fedce3cc5af0618abe520beae617 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 25 Oct 2012 16:24:52 +0000 Subject: [PATCH 08/12] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9007 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/RIGID/fix_rigid.cpp | 64 ++++++++--------------------------------- 1 file changed, 12 insertions(+), 52 deletions(-) diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp index 5c0a7f816d..32a8b7478d 100644 --- a/src/RIGID/fix_rigid.cpp +++ b/src/RIGID/fix_rigid.cpp @@ -695,39 +695,20 @@ void FixRigid::setup(int vflag) tagint *image = atom->image; double **x = atom->x; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - double xy = domain->xy; - double xz = domain->xz; - double yz = domain->yz; + double dx,dy,dz; + double unwrap[3]; for (ibody = 0; ibody < nbody; ibody++) for (i = 0; i < 6; i++) sum[ibody][i] = 0.0; - int xbox,ybox,zbox; - double xunwrap,yunwrap,zunwrap,dx,dy,dz; for (i = 0; i < nlocal; i++) { if (body[i] < 0) continue; ibody = body[i]; - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - - if (triclinic == 0) { - xunwrap = x[i][0] + xbox*xprd; - yunwrap = x[i][1] + ybox*yprd; - zunwrap = x[i][2] + zbox*zprd; - } else { - xunwrap = x[i][0] + xbox*xprd + ybox*xy + zbox*xz; - yunwrap = x[i][1] + ybox*yprd + zbox*yz; - zunwrap = x[i][2] + zbox*zprd; - } - - dx = xunwrap - xcm[ibody][0]; - dy = yunwrap - xcm[ibody][1]; - dz = zunwrap - xcm[ibody][2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - xcm[ibody][0]; + dy = unwrap[1] - xcm[ibody][1]; + dz = unwrap[2] - xcm[ibody][2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; @@ -929,17 +910,9 @@ void FixRigid::final_integrate() double **f = atom->f; int nlocal = atom->nlocal; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - if (triclinic) { - xy = domain->xy; - xz = domain->xz; - yz = domain->yz; - } + double dx,dy,dz; + double unwrap[3]; - int xbox,ybox,zbox; - double xunwrap,yunwrap,zunwrap,dx,dy,dz; for (ibody = 0; ibody < nbody; ibody++) for (i = 0; i < 6; i++) sum[ibody][i] = 0.0; @@ -951,23 +924,10 @@ void FixRigid::final_integrate() sum[ibody][1] += f[i][1]; sum[ibody][2] += f[i][2]; - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - - if (triclinic == 0) { - xunwrap = x[i][0] + xbox*xprd; - yunwrap = x[i][1] + ybox*yprd; - zunwrap = x[i][2] + zbox*zprd; - } else { - xunwrap = x[i][0] + xbox*xprd + ybox*xy + zbox*xz; - yunwrap = x[i][1] + ybox*yprd + zbox*yz; - zunwrap = x[i][2] + zbox*zprd; - } - - dx = xunwrap - xcm[ibody][0]; - dy = yunwrap - xcm[ibody][1]; - dz = zunwrap - xcm[ibody][2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - xcm[ibody][0]; + dy = unwrap[1] - xcm[ibody][1]; + dz = unwrap[2] - xcm[ibody][2]; sum[ibody][3] += dy*f[i][2] - dz*f[i][1]; sum[ibody][4] += dz*f[i][0] - dx*f[i][2]; From e9fa893b1c9735ec4e1d155a5d383648c978ec04 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 25 Oct 2012 16:28:17 +0000 Subject: [PATCH 09/12] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9008 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-MISC/compute_temp_rotate.cpp | 42 +++++++++------------------ src/USER-MISC/fix_addtorque.cpp | 37 +++++++++++------------ 2 files changed, 30 insertions(+), 49 deletions(-) diff --git a/src/USER-MISC/compute_temp_rotate.cpp b/src/USER-MISC/compute_temp_rotate.cpp index b796345404..f76c5fd4b3 100644 --- a/src/USER-MISC/compute_temp_rotate.cpp +++ b/src/USER-MISC/compute_temp_rotate.cpp @@ -89,11 +89,8 @@ double ComputeTempRotate::compute_scalar() { double vthermal[3]; double vcm[3],xcm[3],inertia[3][3],angmom[3],omega[3]; - int xbox,ybox,zbox; double dx,dy,dz; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; + double unwrap[3]; invoked_scalar = update->ntimestep; @@ -123,17 +120,13 @@ double ComputeTempRotate::compute_scalar() for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - xcm[0]; - dy = (x[i][1] + ybox*yprd) - xcm[1]; - dz = (x[i][2] + zbox*zprd) - xcm[2]; + domain->unmap(x[i],image[i],unwrap) + dx = unwrap[0] - xcm[0]; + dy = unwrap[1] - xcm[1]; + dz = unwrap[2] - xcm[2]; vbiasall[i][0] = vcm[0] + dz*omega[1]-dy*omega[2]; vbiasall[i][1] = vcm[1] + dx*omega[2]-dz*omega[0]; vbiasall[i][2] = vcm[2] + dy*omega[0]-dx*omega[1]; - vthermal[0] = v[i][0] - vbiasall[i][0]; vthermal[1] = v[i][1] - vbiasall[i][1]; vthermal[2] = v[i][2] - vbiasall[i][2]; @@ -155,14 +148,10 @@ double ComputeTempRotate::compute_scalar() void ComputeTempRotate::compute_vector() { - int i; double vthermal[3]; double vcm[3],xcm[3],inertia[3][3],angmom[3],omega[3]; - int xbox,ybox,zbox; double dx,dy,dz; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; + double unwrap[3]; invoked_vector = update->ntimestep; @@ -189,25 +178,20 @@ void ComputeTempRotate::compute_vector() } double massone,t[6]; - for (i = 0; i < 6; i++) t[i] = 0.0; + for (int i = 0; i < 6; i++) t[i] = 0.0; - for (i = 0; i < nlocal; i++) + for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - xcm[0]; - dy = (x[i][1] + ybox*yprd) - xcm[1]; - dz = (x[i][2] + zbox*zprd) - xcm[2]; + domain->unmap(x[i],image[i],unwrap) + dx = unwrap[0] - xcm[0]; + dy = unwrap[1] - xcm[1]; + dz = unwrap[2] - xcm[2]; vbiasall[i][0] = vcm[0] + dz*omega[1]-dy*omega[2]; vbiasall[i][1] = vcm[1] + dx*omega[2]-dz*omega[0]; vbiasall[i][2] = vcm[2] + dy*omega[0]-dx*omega[1]; - vthermal[0] = v[i][0] - vbiasall[i][0]; vthermal[1] = v[i][1] - vbiasall[i][1]; vthermal[2] = v[i][2] - vbiasall[i][2]; - if (rmass) massone = rmass[i]; else massone = mass[type[i]]; t[0] += massone * vthermal[0]*vthermal[0]; @@ -219,7 +203,7 @@ void ComputeTempRotate::compute_vector() } MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; + for (int i = 0; i < 6; i++) vector[i] *= force->mvv2e; } /* ---------------------------------------------------------------------- diff --git a/src/USER-MISC/fix_addtorque.cpp b/src/USER-MISC/fix_addtorque.cpp index 7862b2e499..6fd397a4a2 100644 --- a/src/USER-MISC/fix_addtorque.cpp +++ b/src/USER-MISC/fix_addtorque.cpp @@ -109,19 +109,22 @@ void FixAddTorque::init() if (xstr) { xvar = input->variable->find(xstr); - if (xvar < 0) error->all(FLERR,"Variable name for fix addtorque does not exist"); + if (xvar < 0) + error->all(FLERR,"Variable name for fix addtorque does not exist"); if (input->variable->equalstyle(xvar)) xstyle = EQUAL; else error->all(FLERR,"Variable for fix addtorque is invalid style"); } if (ystr) { yvar = input->variable->find(ystr); - if (yvar < 0) error->all(FLERR,"Variable name for fix addtorque does not exist"); + if (yvar < 0) + error->all(FLERR,"Variable name for fix addtorque does not exist"); if (input->variable->equalstyle(yvar)) ystyle = EQUAL; else error->all(FLERR,"Variable for fix addtorque is invalid style"); } if (zstr) { zvar = input->variable->find(zstr); - if (zvar < 0) error->all(FLERR,"Variable name for fix addtorque does not exist"); + if (zvar < 0) + error->all(FLERR,"Variable name for fix addtorque does not exist"); if (input->variable->equalstyle(zvar)) zstyle = EQUAL; else error->all(FLERR,"Variable for fix addtorque is invalid style"); } @@ -168,16 +171,14 @@ void FixAddTorque::post_force(int vflag) int nlocal = atom->nlocal; double mvv2e = force->mvv2e; - int xbox,ybox,zbox; double dx,dy,dz,vx,vy,vz,fx,fy,fz,massone,omegadotr; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - double tcm[3],xcm[3],angmom[3],omega[3],itorque[3],domegadt[3],tlocal[3]; double inertia[3][3]; + double unwrap[3]; + // foriginal[0] = "potential energy" for added force // foriginal[123] = torque on atoms before extra force added + foriginal[0] = foriginal[1] = foriginal[2] = foriginal[3] = 0.0; force_flag = 0; @@ -200,12 +201,10 @@ void FixAddTorque::post_force(int vflag) tlocal[0] = tlocal[1] = tlocal[2] = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - xcm[0]; - dy = (x[i][1] + ybox*yprd) - xcm[1]; - dz = (x[i][2] + zbox*zprd) - xcm[2]; + domain->unmap(x[i],image[i],unwrap) + dx = unwrap[0] - xcm[0]; + dy = unwrap[1] - xcm[1]; + dz = unwrap[2] - xcm[2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; omegadotr = omega[0]*dx+omega[1]*dy+omega[2]*dz; @@ -222,12 +221,10 @@ void FixAddTorque::post_force(int vflag) for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - xcm[0]; - dy = (x[i][1] + ybox*yprd) - xcm[1]; - dz = (x[i][2] + zbox*zprd) - xcm[2]; + domain->unmap(x[i],image[i],unwrap) + dx = unwrap[0] - xcm[0]; + dy = unwrap[1] - xcm[1]; + dz = unwrap[2] - xcm[2]; vx = mvv2e*(dz*omega[1]-dy*omega[2]); vy = mvv2e*(dx*omega[2]-dz*omega[0]); vz = mvv2e*(dy*omega[0]-dx*omega[1]); From dba890a5afae980e586418776eabe08e76b5f5a7 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 25 Oct 2012 16:29:13 +0000 Subject: [PATCH 10/12] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9009 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-MISC/compute_temp_rotate.cpp | 4 ++-- src/USER-MISC/fix_addtorque.cpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/USER-MISC/compute_temp_rotate.cpp b/src/USER-MISC/compute_temp_rotate.cpp index f76c5fd4b3..af5da2ee08 100644 --- a/src/USER-MISC/compute_temp_rotate.cpp +++ b/src/USER-MISC/compute_temp_rotate.cpp @@ -120,7 +120,7 @@ double ComputeTempRotate::compute_scalar() for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - domain->unmap(x[i],image[i],unwrap) + domain->unmap(x[i],image[i],unwrap); dx = unwrap[0] - xcm[0]; dy = unwrap[1] - xcm[1]; dz = unwrap[2] - xcm[2]; @@ -182,7 +182,7 @@ void ComputeTempRotate::compute_vector() for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - domain->unmap(x[i],image[i],unwrap) + domain->unmap(x[i],image[i],unwrap); dx = unwrap[0] - xcm[0]; dy = unwrap[1] - xcm[1]; dz = unwrap[2] - xcm[2]; diff --git a/src/USER-MISC/fix_addtorque.cpp b/src/USER-MISC/fix_addtorque.cpp index 6fd397a4a2..f4dcba3f82 100644 --- a/src/USER-MISC/fix_addtorque.cpp +++ b/src/USER-MISC/fix_addtorque.cpp @@ -201,7 +201,7 @@ void FixAddTorque::post_force(int vflag) tlocal[0] = tlocal[1] = tlocal[2] = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - domain->unmap(x[i],image[i],unwrap) + domain->unmap(x[i],image[i],unwrap); dx = unwrap[0] - xcm[0]; dy = unwrap[1] - xcm[1]; dz = unwrap[2] - xcm[2]; @@ -221,7 +221,7 @@ void FixAddTorque::post_force(int vflag) for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - domain->unmap(x[i],image[i],unwrap) + domain->unmap(x[i],image[i],unwrap); dx = unwrap[0] - xcm[0]; dy = unwrap[1] - xcm[1]; dz = unwrap[2] - xcm[2]; From 7d1c3ed324c08e8cd62203a84ecafc18bd0b3ba9 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 25 Oct 2012 16:33:16 +0000 Subject: [PATCH 11/12] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9010 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/version.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/version.h b/src/version.h index 5d2c676261..57b024e719 100644 --- a/src/version.h +++ b/src/version.h @@ -1 +1 @@ -#define LAMMPS_VERSION "25 Oct 2012" +#define LAMMPS_VERSION "26 Oct 2012" From 8663c9229319c2d6bcf76bb843ccfc8f38a6f72f Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 25 Oct 2012 16:36:45 +0000 Subject: [PATCH 12/12] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9012 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- python/README | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/README b/python/README index 443073f18d..27997582b5 100644 --- a/python/README +++ b/python/README @@ -9,7 +9,7 @@ doc/Section_python.html and in doc/Section_start.html#start_5. Basically you need to follow these steps in the src directory: % make makeshlib # creates Makefile.shlib -% make -f Makefile.shlib g++ # or whatever machine target you wish +% make -f Makefile.shlib g++ # build for whatever machine target you wish % make install-python # may need to do this via sudo You can replace the last step with running the python/install.py