Running NEB calculations through library interface in Fortran

LAMMPS LAMMPS General Discussion
, ,
1

I’m trying to use LAMMPS NEB implementation in my Fortran code. Regular minimization runs fine, but passing the -partition flag necessary for NEB calculations doesn’t seem to not work. I’m using the most recent LAMMPS version, Apr 24, 2023.

If the LAMMPS object is initialized without arguments, i.e. lmp=lammps(), trying to run NEB returns “ERROR: Cannot use NEB with a single replica”, much as expected. Trying to set

args(1)="-partition"
args(2)="11x1" !# example for 11 replicas with 1 core per replica
lmp=lammps(args=args)

fails to read the -partition flag at all, apparently due to how the argv array is constructed in fortran/lammps.f90:lmp_open(). Modifying that a little bit so that argc=size(args)+1 and argv(i+1) = f2c_string(args(i)) passes the flag and the 11x1 correctly, but the next error then is “ERROR: Must use -in switch with multiple partitions”. At this point I ran out of creativity, especially as the contructor “will ignore the -in/-i flag” according to the documentation.

I compiled lammps with make mode=shared mpi, using g++ (Ubuntu 9.4.0-1ubuntu1~20.04.1) 9.4.0. The command to run my fortran code on my laptop is mpirun -np 11 -oversubscribe [executable name]. Any help appreciated.

2

Here is a minimal example for a run using the -partition argument with the fortran library interface:

PROGRAM testlib
  USE LIBLAMMPS                 ! include the LAMMPS library interface
  IMPLICIT NONE
  TYPE(lammps) :: lmp           ! derived type to hold LAMMPS instance
  CHARACTER(LEN=12), PARAMETER :: args(7) = &
      [ CHARACTER(LEN=12) :: 'liblammps', '-log', 'none', '-in', 'none', '-partition', '11x1' ]

  ! create a LAMMPS instance (and initialize MPI)
  lmp = lammps(args)
  ! get and print numerical version code
  PRINT*, 'LAMMPS Version: ', lmp%version()
  ! delete LAMMPS instance (and shutdown MPI)
  CALL lmp%close(.TRUE.)
END PROGRAM testlib

I compile this with:

gfortran -o multi.x -Wall ../fortran/lammps.f90 multi.f90 -L. -llammps -Wl,-rpath,$PWD

And run with:

mpirun -np 11 ./multi.x
1 Like
3

Cheers, this successfully passes the partition flag!

I am now able to execute a NEB calculation from Fortran with equivalent results to the standalone program, but ran into an issue with the program crashing in segmentation fault after NEB finishes. I couldn’t find hints for debugging in the screen files, they all seem to indicate a successfully ended run.

Here is a minimal working example with the aforedescribed behavior. Commands to compile and run are the same as in your example, with the exception that I compile with mpifort instead of gfortran since I’m using the MPI library to call comm rank and size to interpolate the position of jumping atom.

program testlib
  use mpi
  use liblammps                 ! include the LAMMPS library interface
  implicit none
  type(lammps) :: lmp           ! derived type to hold LAMMPS instance
  character(len=12), parameter :: args(7) = &
      [ character(len=12) :: 'liblammps', '-log', 'none', '-in', 'none', '-partition', '11x1' ]
  character(len=512) :: cmds
  integer :: i,myrank,nproc,error
  integer,dimension(18) :: types=1,ids=[(i,i=1,18)]
  !# 3x3x2 fcc {100} slab with nearest neighbour distance equal to 1.0
  double precision,dimension(3*18) :: x=[0.0,0.0,0.0,&
                                         1.0,0.0,0.0,&
                                         2.0,0.0,0.0,&
                                         0.0,1.0,0.0,&
                                         1.0,1.0,0.0,&
                                         2.0,1.0,0.0,&
                                         0.0,2.0,0.0,&
                                         1.0,2.0,0.0,&
                                         2.0,2.0,0.0,&
                                         0.5,0.5,sqrt(0.5),&
                                         1.5,0.5,sqrt(0.5),&
                                         2.5,0.5,sqrt(0.5),&
                                         0.5,1.5,sqrt(0.5),&
                                         1.5,1.5,sqrt(0.5),&
                                         2.5,1.5,sqrt(0.5),&
                                         0.5,2.5,sqrt(0.5),&
                                         1.5,2.5,sqrt(0.5),&
                                         2.5,2.5,sqrt(0.5)]
  !# Initial and final coordinate of a jumping adatom
  double precision,dimension(3) :: xinit=[0.0,0.0,sqrt(2.0)],xfinal=[0.0,1.0,sqrt(2.0)],xjump

  ! create a LAMMPS instance (and initialize MPI)
  lmp = lammps(args)

  cmds="units metal"//new_line("a")//&
       "atom_modify map array sort 0 0.0"//new_line("a")//&
       "region box block 0 3 0 3 0 5"//new_line("a")//&
       "create_box 1 box"//new_line("a")//&
       "pair_style lj/cut 2.5"//new_line("a")//&
       "pair_coeff * * 1 0.890899"//new_line("a")//& !# Lennard-Jonesium with nn distance 0.890899*2^(1/6)=1.0
       "mass 1 1"
  call lmp%commands_string(cmds)

  !# Add slab atoms
  call lmp%create_atoms(id=ids,type=types,x=x)
  !# Freeze the bottom layer
  call lmp%command("group frozen id 1:9")
  call lmp%command("fix freeze frozen setforce 0.0 0.0 0.0")

  !# Add jumping atom; interpolation based on processor rank
  !# to pass replica coordinates without file I/O
  call mpi_comm_rank(MPI_COMM_WORLD,myrank,error)
  call mpi_comm_size(MPI_COMM_WORLD,nproc,error)
  xjump=xinit+myrank*(xfinal-xinit)/(nproc-1)
  call lmp%create_atoms(id=[19],type=[1],x=xjump)

  cmds="fix neb_fix all neb 1.0"//new_line("a")//&
       "min_style quickmin"//new_line("a")//&
       "min_modify dmax 0.01"//new_line("a")//&
       "neb 1e-3 1e-4 1000 1000 100 none"
  call lmp%commands_string(cmds)

  !# Attempt to see if each replica can report its own energy to
  !# get the barrier without file I/O; crashes before this line
  print *,myrank,lmp%get_thermo("pe")

  !# delete LAMMPS instance (and shutdown MPI)
  call lmp%close(.true.)
  
end program testlib
1 Like
4

Please apply this change to the LAMMPS source code, recompile and try again:

  diff --git a/src/REPLICA/neb.cpp b/src/REPLICA/neb.cpp
  index 94b1f860bd..630427b690 100644
  --- a/src/REPLICA/neb.cpp
  +++ b/src/REPLICA/neb.cpp
  @@ -44,7 +44,7 @@ enum { DEFAULT, TERSE, VERBOSE };
   
   /* ---------------------------------------------------------------------- */
   
  -NEB::NEB(LAMMPS *lmp) : Command(lmp), all(nullptr), rdist(nullptr) {}
  +NEB::NEB(LAMMPS *lmp) : Command(lmp), fp(nullptr), all(nullptr), rdist(nullptr) {}
   
   /* ----------------------------------------------------------------------
      internal NEB constructor, called from TAD
1 Like
5

Thanks, works like a charm!

6

It was easy to find thanks to posting a good MWE.
Much appreciated. I wish more people would do that.

1 Like