diff --git a/doc/src/Section_packages.txt b/doc/src/Section_packages.txt
index 18030162cf3f9666c73ef9b7862d6e39895ac601..76f88b8ab81967f53e9dd96cb6402bb7dbcbc5e2 100644
--- a/doc/src/Section_packages.txt
+++ b/doc/src/Section_packages.txt
@@ -2062,8 +2062,9 @@ to plain C++. In contrast to the MEAM package, no library
 needs to be compiled and the pair style can be instantiated
 multiple times.
 
-[Author:]  Sebastian Huetter, (Otto-von-Guericke University Magdeburg)
-based on the work of Greg Wagner (Northwestern U) while at Sandia.
+[Author:] Sebastian Huetter, (Otto-von-Guericke University Magdeburg)
+based on the Fortran version of Greg Wagner (Northwestern U) while at
+Sandia.
 
 [Install or un-install:]
   
diff --git a/examples/COUPLE/README b/examples/COUPLE/README
index 7a62536239d3d8af31265f2790f3fc2e0a112da4..c8c9e0e31b0e350a4ef9af6478bf973b11cb9493 100644
--- a/examples/COUPLE/README
+++ b/examples/COUPLE/README
@@ -41,5 +41,8 @@ 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
+fortran3            wrapper written by Nir Goldman (LLNL), as an
+                      extension to fortran2, used for calling LAMMPS
+                      from Fortran DFTB+ code
 
-Each sub-directory has its own README.
+Each sub-directory has its own README with more details.
diff --git a/examples/COUPLE/fortran3/LAMMPS-wrapper.cpp b/examples/COUPLE/fortran3/LAMMPS-wrapper.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6e8bbec5ae292198153079289568c16408103ff7
--- /dev/null
+++ b/examples/COUPLE/fortran3/LAMMPS-wrapper.cpp
@@ -0,0 +1,236 @@
+/* -----------------------------------------------------------------------
+    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 <karlh@ugcs.caltech.edu>
+                          University of Tennessee, Knoxville (USA), 2012
+------------------------------------------------------------------------- */
+
+/* This is set of "wrapper" functions to assist LAMMPS.F90, which itself
+   provides a (I hope) robust Fortran interface to library.cpp and
+   library.h.  All functions herein COULD be added to library.cpp instead of
+   including this as a separate file. See the README for instructions. */
+
+#include <mpi.h>
+#include "LAMMPS-wrapper.h"
+#include <library.h>
+#include <lammps.h>
+#include <atom.h>
+#include <fix.h>
+#include <compute.h>
+#include <modify.h>
+#include <error.h>
+#include <cstdlib>
+
+using namespace LAMMPS_NS;
+
+void lammps_open_fortran_wrapper (int argc, char **argv,
+      MPI_Fint communicator, void **ptr)
+{
+   MPI_Comm C_communicator = MPI_Comm_f2c (communicator);
+   lammps_open (argc, argv, C_communicator, ptr);
+}
+
+int lammps_get_ntypes (void *ptr)
+{
+  class LAMMPS *lmp = (class LAMMPS *) ptr;
+  int ntypes = lmp->atom->ntypes;
+  return ntypes;
+}
+
+void lammps_error_all (void *ptr, const char *file, int line, const char *str)
+{
+   class LAMMPS *lmp = (class LAMMPS *) ptr;
+   lmp->error->all (file, line, str);
+}
+
+int lammps_extract_compute_vectorsize (void *ptr, char *id, int style)
+{
+   class LAMMPS *lmp = (class LAMMPS *) ptr;
+   int icompute = lmp->modify->find_compute(id);
+   if ( icompute < 0 ) return 0;
+   class Compute *compute = lmp->modify->compute[icompute];
+
+   if ( style == 0 )
+   {
+      if ( !compute->vector_flag )
+         return 0;
+      else
+         return compute->size_vector;
+   }
+   else if ( style == 1 )
+   {
+      return lammps_get_natoms (ptr);
+   }
+   else if ( style == 2 )
+   {
+      if ( !compute->local_flag )
+         return 0;
+      else
+         return compute->size_local_rows;
+   }
+   else
+      return 0;
+}
+
+void lammps_extract_compute_arraysize (void *ptr, char *id, int style,
+      int *nrows, int *ncols)
+{
+   class LAMMPS *lmp = (class LAMMPS *) ptr;
+   int icompute = lmp->modify->find_compute(id);
+   if ( icompute < 0 )
+   {
+      *nrows = 0;
+      *ncols = 0;
+   }
+   class Compute *compute = lmp->modify->compute[icompute];
+
+   if ( style == 0 )
+   {
+      if ( !compute->array_flag )
+      {
+         *nrows = 0;
+         *ncols = 0;
+      }
+      else
+      {
+         *nrows = compute->size_array_rows;
+         *ncols = compute->size_array_cols;
+      }
+   }
+   else if ( style == 1 )
+   {
+      if ( !compute->peratom_flag )
+      {
+         *nrows = 0;
+         *ncols = 0;
+      }
+      else
+      {
+         *nrows = lammps_get_natoms (ptr);
+         *ncols = compute->size_peratom_cols;
+      }
+   }
+   else if ( style == 2 )
+   {
+      if ( !compute->local_flag )
+      {
+         *nrows = 0;
+         *ncols = 0;
+      }
+      else
+      {
+         *nrows = compute->size_local_rows;
+         *ncols = compute->size_local_cols;
+      }
+   }
+   else
+   {
+      *nrows = 0;
+      *ncols = 0;
+   }
+
+   return;
+}
+
+int lammps_extract_fix_vectorsize (void *ptr, char *id, int style)
+{
+   class LAMMPS *lmp = (class LAMMPS *) ptr;
+   int ifix = lmp->modify->find_fix(id);
+   if ( ifix < 0 ) return 0;
+   class Fix *fix = lmp->modify->fix[ifix];
+
+   if ( style == 0 )
+   {
+      if ( !fix->vector_flag )
+         return 0;
+      else
+         return fix->size_vector;
+   }
+   else if ( style == 1 )
+   {
+      return lammps_get_natoms (ptr);
+   }
+   else if ( style == 2 )
+   {
+      if ( !fix->local_flag )
+         return 0;
+      else
+         return fix->size_local_rows;
+   }
+   else
+      return 0;
+}
+
+void lammps_extract_fix_arraysize (void *ptr, char *id, int style,
+      int *nrows, int *ncols)
+{
+   class LAMMPS *lmp = (class LAMMPS *) ptr;
+   int ifix = lmp->modify->find_fix(id);
+   if ( ifix < 0 )
+   {
+      *nrows = 0;
+      *ncols = 0;
+   }
+   class Fix *fix = lmp->modify->fix[ifix];
+
+   if ( style == 0 )
+   {
+      if ( !fix->array_flag )
+      {
+         *nrows = 0;
+         *ncols = 0;
+      }
+      else
+      {
+         *nrows = fix->size_array_rows;
+         *ncols = fix->size_array_cols;
+      }
+   }
+   else if ( style == 1 )
+   {
+      if ( !fix->peratom_flag )
+      {
+         *nrows = 0;
+         *ncols = 0;
+      }
+      else
+      {
+         *nrows = lammps_get_natoms (ptr);
+         *ncols = fix->size_peratom_cols;
+      }
+   }
+   else if ( style == 2 )
+   {
+      if ( !fix->local_flag )
+      {
+         *nrows = 0;
+         *ncols = 0;
+      }
+      else
+      {
+         *nrows = fix->size_local_rows;
+         *ncols = fix->size_local_cols;
+      }
+   }
+   else
+   {
+      *nrows = 0;
+      *ncols = 0;
+   }
+
+   return;
+
+}
+
+/* vim: set ts=3 sts=3 expandtab: */
diff --git a/examples/COUPLE/fortran3/LAMMPS-wrapper.h b/examples/COUPLE/fortran3/LAMMPS-wrapper.h
new file mode 100644
index 0000000000000000000000000000000000000000..7d94362436f82f1e8a9bdc6e51d7bb71308c89bd
--- /dev/null
+++ b/examples/COUPLE/fortran3/LAMMPS-wrapper.h
@@ -0,0 +1,40 @@
+/* -----------------------------------------------------------------------
+    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 <karlh@ugcs.caltech.edu>
+                          University of Tennessee, Knoxville (USA), 2012
+------------------------------------------------------------------------- */
+
+/* This is set of "wrapper" functions to assist LAMMPS.F90, which itself
+   provides a (I hope) robust Fortran interface to library.cpp and
+   library.h.  All prototypes herein COULD be added to library.h instead of
+   including this as a separate file. See the README for instructions. */
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/* Prototypes for auxiliary functions */
+void lammps_open_fortran_wrapper (int, char**, MPI_Fint, void**);
+int lammps_get_ntypes (void*);
+int lammps_extract_compute_vectorsize (void*, char*, int);
+void lammps_extract_compute_arraysize (void*, char*, int, int*, int*);
+int lammps_extract_fix_vectorsize (void*, char*, int);
+void lammps_extract_fix_arraysize (void*, char*, int, int*, int*);
+void lammps_error_all (void*, const char*, int, const char*);
+
+#ifdef __cplusplus
+}
+#endif
+
+/* vim: set ts=3 sts=3 expandtab: */
diff --git a/examples/COUPLE/fortran3/LAMMPS-wrapper2.cpp b/examples/COUPLE/fortran3/LAMMPS-wrapper2.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f245c44d79c26fd13266b8b5e301c427b114cb73
--- /dev/null
+++ b/examples/COUPLE/fortran3/LAMMPS-wrapper2.cpp
@@ -0,0 +1,57 @@
+/* -----------------------------------------------------------------------
+    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 <karlh@ugcs.caltech.edu>
+                          University of Tennessee, Knoxville (USA), 2012
+------------------------------------------------------------------------- */
+
+/* This is set of "wrapper" functions to assist LAMMPS.F90, which itself
+   provides a (I hope) robust Fortran interface to library.cpp and
+   library.h.  All functions herein COULD be added to library.cpp instead of
+   including this as a separate file. See the README for instructions. */
+
+#include <mpi.h>
+#include "LAMMPS-wrapper2.h"
+#include <library.h>
+#include <lammps.h>
+#include <atom.h>
+#include <input.h>
+#include <modify.h>
+#include <fix.h>
+#include <fix_external.h>
+#include <compute.h>
+#include <modify.h>
+#include <error.h>
+#include <cstdlib>
+
+using namespace LAMMPS_NS;
+
+extern "C" void f_callback(void *, bigint, int, tagint *, double **, double **);
+
+void lammps_set_callback (void *ptr) {
+  class LAMMPS *lmp = (class LAMMPS *) ptr;
+  int ifix = lmp->modify->find_fix_by_style("external");
+  FixExternal *fix = (FixExternal *) lmp->modify->fix[ifix];
+  fix->set_callback(f_callback, ptr);
+  return;
+}
+
+void lammps_set_user_energy (void *ptr, double energy) {
+  class LAMMPS *lmp = (class LAMMPS *) ptr;
+  int ifix = lmp->modify->find_fix_by_style("external");
+  FixExternal *fix = (FixExternal *) lmp->modify->fix[ifix];
+  fix->set_energy(energy);
+  return;
+}
+
diff --git a/examples/COUPLE/fortran3/LAMMPS-wrapper2.h b/examples/COUPLE/fortran3/LAMMPS-wrapper2.h
new file mode 100644
index 0000000000000000000000000000000000000000..794006e3afac4919d33cf1d0080b1d4a0a22dc49
--- /dev/null
+++ b/examples/COUPLE/fortran3/LAMMPS-wrapper2.h
@@ -0,0 +1,34 @@
+/* -----------------------------------------------------------------------
+    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:  Nir Goldman, ngoldman@llnl.gov, Oct. 19th, 2016
+------------------------------------------------------------------------- */
+
+/* This is set of "wrapper" functions to assist LAMMPS.F90, which itself
+   provides a (I hope) robust Fortran interface to library.cpp and
+   library.h.  All prototypes herein COULD be added to library.h instead of
+   including this as a separate file. See the README for instructions. */
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/* Prototypes for auxiliary functions */
+void lammps_set_callback (void *); 
+void lammps_set_user_energy (void*, double); 
+
+#ifdef __cplusplus
+}
+#endif
+
+/* vim: set ts=3 sts=3 expandtab: */
diff --git a/examples/COUPLE/fortran3/LAMMPS.F90 b/examples/COUPLE/fortran3/LAMMPS.F90
new file mode 100644
index 0000000000000000000000000000000000000000..eb5b7f825becfb54fd80c2400c3f4ef52fd94b82
--- /dev/null
+++ b/examples/COUPLE/fortran3/LAMMPS.F90
@@ -0,0 +1,956 @@
+!! -----------------------------------------------------------------------
+!   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 <karlh@ugcs.caltech.edu>
+!                         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, lammps_set_callback, lammps_set_user_energy
+   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_set_callback (ptr) &
+      bind (C, name='lammps_set_callback')
+        import :: C_ptr
+        type (C_ptr), value :: ptr
+      end subroutine lammps_set_callback
+
+      subroutine lammps_set_user_energy (ptr, energy) &
+      bind (C, name='lammps_set_user_energy')
+        import :: C_ptr, C_double 
+        type (C_ptr), value :: ptr
+        real(C_double), value :: energy
+      end subroutine lammps_set_user_energy 
+
+      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
+         atom => 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' .and. name /= 'fexternal' ) 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)
+      Cptr = lammps_extract_compute_Cptr (ptr, id, style, type)
+      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/fortran3/README b/examples/COUPLE/fortran3/README
new file mode 100644
index 0000000000000000000000000000000000000000..9effa35ec4e0040bcfca794be6757aad05e28fd2
--- /dev/null
+++ b/examples/COUPLE/fortran3/README
@@ -0,0 +1,33 @@
+This directory has an example of using a callback function to obtain
+forces from a fortran code for a LAMMPS simulation.  The reader should
+refer to the README file in COUPLE/fortran2 before proceeding. Here,
+the LAMMPS.F90 file has been modified slightly and additional files
+named LAMMPS-wrapper2.h and LAMMPS-wrapper2.cpp have been included in
+order to supply wrapper functions to set the LAMMPS callback function
+and total energy.
+
+In this example, the callback function is set to run the
+semi-empirical quantum code DFTB+ in serial and then read in the total
+energy, forces, and stress tensor from file.  In this case, nlocal =
+the total number of atoms in the system, so particle positions can be
+read from the pos array directly, and DFTB+ forces can simply be
+included via the fext array. The user should take care in the case of
+a parallel calculation, where LAMMPS can assign different particules
+to each processor. For example, the user should use functions such as
+lammps_gather_atoms() and lammps_scatter_atoms() in the case where the
+fortran force calculating code requires the positions of all atoms,
+etc.
+
+A few more important notes: 
+
+-The stress tensor from DFTB+ is passed in to LAMMPS via pointer.
+-Calling the subroutine lammps_set_callback() is required in order to set
+       a pointer to the callback function in LAMMPS.
+-The subroutine lammps_set_user_energy() passes in the potential energy
+       from DFTB+ to LAMMPS.
+
+This example was created by Nir Goldman, whom you can contact with
+questions:
+
+Nir Goldman, LLNL
+ngoldman@llnl.gov
diff --git a/examples/COUPLE/fortran3/data.diamond b/examples/COUPLE/fortran3/data.diamond
new file mode 100644
index 0000000000000000000000000000000000000000..b3dd599cf4755ec2bef6eed4094dd38b9e34e200
--- /dev/null
+++ b/examples/COUPLE/fortran3/data.diamond
@@ -0,0 +1,148 @@
+# Position data file
+
+64 atoms
+1 atom types
+
+0 7.134 xlo xhi
+0 7.134 ylo yhi
+0 7.134 zlo zhi
+
+0.00000000 0.00000000 0.00000000 xy xz yz
+
+Masses
+
+1 12.010000
+
+Atoms
+
+     1 1 0     0 0 0
+     2 1 0     0.89175 0.89175 0.89175
+     3 1 0     1.7835 1.7835 0
+     4 1 0     2.67525 2.67525 0.89175
+     5 1 0     0 1.7835 1.7835
+     6 1 0     0.89175 2.67525 2.67525
+     7 1 0     1.7835 0 1.7835
+     8 1 0     2.67525 0.89175 2.67525
+     9 1 0     0 0 3.567
+     10 1 0     0.89175 0.89175 4.45875
+     11 1 0     1.7835 1.7835 3.567
+     12 1 0     2.67525 2.67525 4.45875
+     13 1 0     0 1.7835 5.3505
+     14 1 0     0.89175 2.67525 6.24225
+     15 1 0     1.7835 0 5.3505
+     16 1 0     2.67525 0.89175 6.24225
+     17 1 0     0 3.567 0
+     18 1 0     0.89175 4.45875 0.89175
+     19 1 0     1.7835 5.3505 0
+     20 1 0     2.67525 6.24225 0.89175
+     21 1 0     0 5.3505 1.7835
+     22 1 0     0.89175 6.24225 2.67525
+     23 1 0     1.7835 3.567 1.7835
+     24 1 0     2.67525 4.45875 2.67525
+     25 1 0     0 3.567 3.567
+     26 1 0     0.89175 4.45875 4.45875
+     27 1 0     1.7835 5.3505 3.567
+     28 1 0     2.67525 6.24225 4.45875
+     29 1 0     0 5.3505 5.3505
+     30 1 0     0.89175 6.24225 6.24225
+     31 1 0     1.7835 3.567 5.3505
+     32 1 0     2.67525 4.45875 6.24225
+     33 1 0     3.567 0 0
+     34 1 0     4.45875 0.89175 0.89175
+     35 1 0     5.3505 1.7835 0
+     36 1 0     6.24225 2.67525 0.89175
+     37 1 0     3.567 1.7835 1.7835
+     38 1 0     4.45875 2.67525 2.67525
+     39 1 0     5.3505 0 1.7835
+     40 1 0     6.24225 0.89175 2.67525
+     41 1 0     3.567 0 3.567
+     42 1 0     4.45875 0.89175 4.45875
+     43 1 0     5.3505 1.7835 3.567
+     44 1 0     6.24225 2.67525 4.45875
+     45 1 0     3.567 1.7835 5.3505
+     46 1 0     4.45875 2.67525 6.24225
+     47 1 0     5.3505 0 5.3505
+     48 1 0     6.24225 0.89175 6.24225
+     49 1 0     3.567 3.567 0
+     50 1 0     4.45875 4.45875 0.89175
+     51 1 0     5.3505 5.3505 0
+     52 1 0     6.24225 6.24225 0.89175
+     53 1 0     3.567 5.3505 1.7835
+     54 1 0     4.45875 6.24225 2.67525
+     55 1 0     5.3505 3.567 1.7835
+     56 1 0     6.24225 4.45875 2.67525
+     57 1 0     3.567 3.567 3.567
+     58 1 0     4.45875 4.45875 4.45875
+     59 1 0     5.3505 5.3505 3.567
+     60 1 0     6.24225 6.24225 4.45875
+     61 1 0     3.567 5.3505 5.3505
+     62 1 0     4.45875 6.24225 6.24225
+     63 1 0     5.3505 3.567 5.3505
+     64 1 0     6.24225 4.45875 6.24225
+
+Velocities
+
+1 -0.00733742 -0.0040297 -0.00315229 
+2 -0.00788609 -0.00567535 -0.00199152 
+3 -0.00239042 0.00710139 -0.00335049 
+4 0.00678551 0.0019976 0.00219289 
+5 0.00413717 0.00275709 0.000937637 
+6 -0.00126313 0.00485636 0.00727862 
+7 0.00337547 -0.00234623 -0.000922223 
+8 -0.00792183 -0.00509186 -0.00104168 
+9 0.00414091 0.00390285 0.000845961 
+10 -0.000284543 0.0010771 -0.00458404 
+11 -0.00394968 -0.00446363 -0.00361688 
+12 0.00067088 -0.00655175 -0.00752464 
+13 0.00306632 -0.00245545 -0.00183867 
+14 -0.0082145 -0.00564127 0.000281191 
+15 0.00504454 0.0045835 0.000495763 
+16 0.0035767 0.00320441 -0.00486426 
+17 0.00420597 0.00262005 -0.0049459 
+18 0.00440579 -1.76783e-05 0.00449311 
+19 -0.00406463 0.00613304 0.00285599 
+20 0.00171215 -0.00517887 0.00124326 
+21 0.0011118 0.00334129 -0.0015222 
+22 -0.00838394 -0.00112906 -0.00353379 
+23 -0.00578527 -0.00415501 0.00297043 
+24 -0.00211466 0.000964108 -0.00716523 
+25 -0.000204107 -0.00380986 0.00681648 
+26 0.00677838 0.00540935 0.0044354 
+27 -0.00266809 -0.00358382 -0.00241889 
+28 -0.0003973 0.00236566 0.00558871 
+29 0.000754103 0.00457797 0.000105531 
+30 -0.00246049 0.00110428 0.00511088 
+31 0.00248891 0.00623314 0.00461597 
+32 -0.00509423 0.000570503 0.00720856 
+33 -0.00244427 -0.00374384 0.00618767 
+34 -0.000360752 -8.10558e-05 0.00314052 
+35 0.00435313 -0.00630587 -0.0070309 
+36 0.00651087 -0.00389833 3.72525e-05 
+37 0.00631828 -0.00316064 0.00231522 
+38 -0.00579624 -0.00345068 -0.000277486 
+39 0.00483974 0.000715028 0.000206355 
+40 -0.00388164 -0.00189242 -0.00554862 
+41 0.00398115 0.00152915 0.00756919 
+42 -0.000552263 0.00352025 -0.000246143 
+43 -0.00800284 0.00555703 0.00425716 
+44 -0.00734405 -0.00752512 0.00667173 
+45 -0.00545636 0.00421035 0.00399552 
+46 0.00480246 0.00621147 -0.00492715 
+47 -0.00424168 0.00621818 -9.37733e-05 
+48 -0.00649561 0.00612908 -0.0020753 
+49 -0.0075007 -0.00384737 -0.00687913 
+50 -0.00203903 -0.00764372 0.0023883 
+51 0.00442642 0.00744072 -0.0049344 
+52 -0.00280486 -0.00509128 -0.00678045 
+53 0.00679491 0.00583493 0.00333875 
+54 0.00574665 -0.00521074 0.00523475 
+55 0.00305618 -0.00320094 0.00341297 
+56 0.004304 0.000615544 -0.00668787 
+57 0.00564532 0.00327373 0.00388611 
+58 0.000676899 0.00210326 0.00495295 
+59 0.000160781 -0.00744313 -0.00279828 
+60 0.00623521 0.00371301 0.00178015 
+61 0.00520759 0.000642669 0.00207913 
+62 0.00398042 0.0046438 -0.00359978 
+63 -0.00478071 -0.00304932 -0.00765125 
+64 0.00282671 -0.00548392 -0.00692691 
diff --git a/examples/COUPLE/fortran3/in.simple b/examples/COUPLE/fortran3/in.simple
new file mode 100644
index 0000000000000000000000000000000000000000..894a490cf8b032b71bbdab9d073a9dee6a40c672
--- /dev/null
+++ b/examples/COUPLE/fortran3/in.simple
@@ -0,0 +1,16 @@
+units real
+atom_style	charge
+atom_modify map array
+atom_modify sort 0 0.0
+read_data data.diamond
+neighbor        1.0 bin
+neigh_modify    delay 0 every 5 check no
+fix 1 all nve
+fix 2 all external pf/callback 1 1
+
+fix_modify 2 energy yes 
+thermo_style custom step temp etotal ke pe lx ly lz pxx pyy pzz press 
+
+thermo          1
+timestep        0.5
+
diff --git a/examples/COUPLE/fortran3/makefile b/examples/COUPLE/fortran3/makefile
new file mode 100644
index 0000000000000000000000000000000000000000..86dea30850676c3a58fee935ca6da612e7c2476d
--- /dev/null
+++ b/examples/COUPLE/fortran3/makefile
@@ -0,0 +1,45 @@
+SHELL = /bin/sh
+
+# Path to LAMMPS extraction directory
+LAMMPS_ROOT = ../../..
+LAMMPS_SRC = $(LAMMPS_ROOT)/src
+
+# Uncomment the line below if using the MPI stubs library
+MPI_STUBS = #-I$(LAMMPS_SRC)/STUBS
+
+FC = mpif90 # replace with your Fortran compiler
+CXX = mpicc # replace with your C++ compiler
+
+# Flags for Fortran compiler, C++ compiler, and C preprocessor, respectively
+FFLAGS = -O2 -fPIC 
+CXXFLAGS = -O2 -fPIC
+CPPFLAGS = -DOMPI_SKIP_MPICXX=1 -DMPICH_SKIP_MPICXX
+
+all : liblammps_fortran.a liblammps_fortran.so simpleF.x
+
+liblammps_fortran.so : LAMMPS.o LAMMPS-wrapper.o LAMMPS-wrapper2.o
+	$(FC) $(FFLAGS) -shared -o $@ $^
+
+simpleF.x: simple.o LAMMPS.o LAMMPS-wrapper.o LAMMPS-wrapper2.o
+	$(FC) $(FFLAGS) simple.o -o simpleF.x liblammps_fortran.a $(LAMMPS_SRC)/liblammps_mvapich.a -lstdc++ /usr/local/tools/fftw/lib/libfftw.a
+
+liblammps_fortran.a : LAMMPS.o LAMMPS-wrapper.o LAMMPS-wrapper2.o
+	$(AR) rs $@ $^
+
+LAMMPS.o lammps.mod : LAMMPS.F90
+	$(FC) $(CPPFLAGS) $(FFLAGS) -c $<
+
+simple.o : simple.f90
+	$(FC) $(FFLAGS) -c $<
+
+LAMMPS-wrapper.o : LAMMPS-wrapper.cpp LAMMPS-wrapper.h
+	$(CXX) $(CPPFLAGS) $(CXXFLAGS) -c $< -I$(LAMMPS_SRC) $(MPI_STUBS)
+
+LAMMPS-wrapper2.o : LAMMPS-wrapper2.cpp LAMMPS-wrapper2.h
+	$(CXX) $(CPPFLAGS) $(CXXFLAGS) -c $< -I$(LAMMPS_SRC) $(MPI_STUBS)
+
+clean :
+	$(RM) *.o *.mod liblammps_fortran.a liblammps_fortran.so
+
+dist :
+	tar -czvf fortran-interface-callback.tar.gz LAMMPS-wrapper.h LAMMPS-wrapper.cpp LAMMPS-wrapper2.h LAMMPS-wrapper2.cpp LAMMPS.F90 makefile README simple.f90
diff --git a/examples/COUPLE/fortran3/simple.f90 b/examples/COUPLE/fortran3/simple.f90
new file mode 100644
index 0000000000000000000000000000000000000000..40f8bf8b86dedb840becd448ad367df30a56a5b0
--- /dev/null
+++ b/examples/COUPLE/fortran3/simple.f90
@@ -0,0 +1,114 @@
+   module callback
+     implicit none
+     contains
+       subroutine fortran_callback(lmp, timestep, nlocal, ids, c_pos, c_fext) & 
+     & bind(C, name='f_callback')
+       use, intrinsic :: ISO_C_binding
+       use LAMMPS
+       implicit none
+       type (C_ptr), value :: lmp
+       integer(C_int64_t), intent(in), value :: timestep
+       integer(C_int), intent(in), value :: nlocal
+       real (C_double), dimension(:,:), pointer :: x
+       type(c_ptr) :: c_pos, c_fext, c_ids
+       double precision, pointer :: fext(:,:), pos(:,:)
+       integer, intent(in) :: ids(nlocal)
+       real (C_double), dimension(:), pointer :: virial => NULL()
+       real (C_double) :: etot
+       real(C_double), pointer :: ts_lmp
+       double precision :: stress(3,3), ts_dftb
+       integer :: natom , i
+       real (C_double), parameter :: econv = 627.4947284155114 ! converts from Ha to
+       double precision, parameter :: fconv = 1185.793095983065 ! converts from Ha/bohr to
+       double precision, parameter :: autoatm = 2.9037166638E8
+       double precision lx, ly, lz
+       real (C_double), pointer :: boxxlo, boxxhi
+       real (C_double), pointer :: boxylo, boxyhi
+       real (C_double), pointer :: boxzlo, boxzhi
+       double precision, parameter :: nktv2p = 68568.4149999999935972
+       double precision :: volume
+       type (C_ptr) :: Cptr
+       type (C_ptr), pointer, dimension(:) :: Catom
+
+       call c_f_pointer(c_pos, pos, [3,nlocal])
+       call c_f_pointer(c_fext, fext, [3,nlocal])
+       call lammps_extract_global(boxxlo, lmp, 'boxxlo')
+       call lammps_extract_global(boxxhi, lmp, 'boxxhi')
+       call lammps_extract_global(boxylo, lmp, 'boxylo')
+       call lammps_extract_global(boxyhi, lmp, 'boxyhi')
+       call lammps_extract_global(boxzlo, lmp, 'boxzlo')
+       call lammps_extract_global(boxzhi, lmp, 'boxzhi')
+       lx = boxxhi - boxxlo
+       ly = boxyhi - boxylo
+       lz = boxzhi - boxzlo
+       volume = lx*ly*lz
+       open (unit = 10, status = 'replace', action = 'write', file='lammps.gen')
+       write(10,*)nlocal,"S"
+       write(10,*) "C"
+       do i = 1, nlocal
+         write(10,'(2I,3F15.6)')i,1,pos(:,ids(i))
+       enddo
+       write(10,*)"0.0 0.0 0.0"
+       write(10,*)lx,0,0
+       write(10,*)0,ly,0
+       write(10,*)0,0,lz
+       close(10)
+       call system("./dftb+ > dftb.out")
+       open (unit = 10, status = 'old', file = 'results.out')
+       read(10,*)etot
+       read(10,*)ts_dftb
+       do i = 1, 3
+         read(10,*)stress(i,:)
+       enddo
+       stress (:,:) = stress(:,:)*autoatm
+       etot = etot*econv
+       call lammps_extract_global(ts_lmp, lmp, 'TS_dftb')
+       ts_lmp = ts_dftb
+       do i = 1, nlocal
+         read(10,*)fext(:,ids(i))
+         fext(:,ids(i)) = fext(:,ids(i))*fconv
+       enddo
+       close(10)
+       call lammps_set_user_energy (lmp, etot)
+       call lammps_extract_atom (virial, lmp, 'virial')
+       if (.not. associated(virial)) then
+         print*,'virial pointer not associated.'
+         STOP
+       endif
+       virial(1) = stress(1,1)/(nktv2p/volume)
+       virial(2) = stress(2,2)/(nktv2p/volume)
+       virial(3) = stress(3,3)/(nktv2p/volume)
+       virial(4) = stress(1,2)/(nktv2p/volume)
+       virial(5) = stress(1,3)/(nktv2p/volume)
+       virial(6) = stress(2,3)/(nktv2p/volume)
+
+       end  subroutine
+    end module callback
+
+
+program simple_fortran_callback
+
+   use MPI
+   use LAMMPS
+   use callback
+   use, intrinsic :: ISO_C_binding, only : C_double, C_ptr, C_int, C_FUNPTR
+   implicit none
+   type (C_ptr) :: lmp
+   integer :: error, narg, me, nprocs
+
+   call MPI_Init (error)
+   call MPI_Comm_rank (MPI_COMM_WORLD, me, error)
+   call MPI_Comm_size (MPI_COMM_WORLD, nprocs, error)
+
+   call lammps_open_no_mpi ('lmp -log log.simple', lmp)
+   call lammps_file (lmp, 'in.simple')
+   call lammps_set_callback(lmp)
+
+   call lammps_command (lmp, 'run 10')
+   call lammps_close (lmp)
+   call MPI_Finalize (error)
+
+
+end program simple_fortran_callback
+
+
diff --git a/src/modify.cpp b/src/modify.cpp
index 01de6b59284676069de3980362accd8aa35badd6..89c4a43001cc24d136a796eb2be4038203302f74 100644
--- a/src/modify.cpp
+++ b/src/modify.cpp
@@ -996,7 +996,6 @@ int Modify::check_package(const char *package_fix_name)
   return 1;
 }
 
-
 /* ----------------------------------------------------------------------
    check if the group indicated by groupbit overlaps with any
    currently existing rigid fixes. return 1 in this case otherwise 0