diff --git a/doc/src/variable.rst b/doc/src/variable.rst
index 6a6cdf9a6d..af324c180f 100644
--- a/doc/src/variable.rst
+++ b/doc/src/variable.rst
@@ -65,8 +65,8 @@ Syntax
                            bound(group,dir,region), gyration(group,region), ke(group,reigon),
                            angmom(group,dim,region), torque(group,dim,region),
                            inertia(group,dimdim,region), omega(group,dim,region)
-         special functions = sum(x), min(x), max(x), ave(x), trap(x), slope(x), gmask(x), rmask(x), grmask(x,y), next(x)
-         feature functions = is_active(category,feature,exact), is_defined(category,id,exact)
+         special functions = sum(x), min(x), max(x), ave(x), trap(x), slope(x), gmask(x), rmask(x), grmask(x,y), next(x), is_file(name)
+         feature functions = is_available(category,feature), is_active(category,feature), is_defined(category,id)
          atom value = id[i], mass[i], type[i], mol[i], x[i], y[i], z[i], vx[i], vy[i], vz[i], fx[i], fy[i], fz[i], q[i]
          atom vector = id, mass, type, mol, x, y, z, vx, vy, vz, fx, fy, fz, q
          compute references = c_ID, c_ID[i], c_ID[i][j], C_ID, C_ID[i]
@@ -821,6 +821,10 @@ Special Functions
 Special functions take specific kinds of arguments, meaning their
 arguments cannot be formulas themselves.
 
+The is_file(x) function is a test whether 'x' is a (readable) file
+and returns 1 in this case, otherwise it returns 0.  For that 'x'
+is taken as a literal string and must not have any blanks in it.
+
 The sum(x), min(x), max(x), ave(x), trap(x), and slope(x) functions
 each take 1 argument which is of the form "c_ID" or "c_ID[N]" or
 "f_ID" or "f_ID[N]" or "v_name".  The first two are computes and the
diff --git a/src/variable.cpp b/src/variable.cpp
index 3bf88db55e..1a4feb3573 100644
--- a/src/variable.cpp
+++ b/src/variable.cpp
@@ -65,7 +65,7 @@ enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,MODULO,UNARY,
      SQRT,EXP,LN,LOG,ABS,SIN,COS,TAN,ASIN,ACOS,ATAN,ATAN2,
      RANDOM,NORMAL,CEIL,FLOOR,ROUND,RAMP,STAGGER,LOGFREQ,LOGFREQ2,
      LOGFREQ3,STRIDE,STRIDE2,VDISPLACE,SWIGGLE,CWIGGLE,GMASK,RMASK,
-     GRMASK,IS_ACTIVE,IS_DEFINED,IS_AVAILABLE,
+     GRMASK,IS_ACTIVE,IS_DEFINED,IS_AVAILABLE,IS_FILE,
      VALUE,ATOMARRAY,TYPEARRAY,INTARRAY,BIGINTARRAY,VECTORARRAY};
 
 // customize by adding a special function
@@ -4079,7 +4079,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
       strcmp(word,"gmask") && strcmp(word,"rmask") &&
       strcmp(word,"grmask") && strcmp(word,"next") &&
       strcmp(word,"is_active") && strcmp(word,"is_defined") &&
-      strcmp(word,"is_available"))
+      strcmp(word,"is_available") && strcmp(word,"is_file"))
     return 0;
 
   // parse contents for comma-separated args
@@ -4488,6 +4488,26 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
 
     // save value in tree or on argstack
 
+    if (tree) {
+      Tree *newtree = new Tree();
+      newtree->type = VALUE;
+      newtree->value = value;
+      newtree->first = newtree->second = nullptr;
+      newtree->nextra = 0;
+      treestack[ntreestack++] = newtree;
+    } else argstack[nargstack++] = value;
+
+  } else if (strcmp(word,"is_file") == 0) {
+    if (narg != 1)
+      print_var_error(FLERR,"Invalid is_file() function in "
+                      "variable formula",ivar);
+
+    FILE *fp = fopen(args[0],"r");
+    value = (fp == nullptr) ? 0.0 : 1.0;
+    if (fp) fclose(fp);
+
+    // save value in tree or on argstack
+
     if (tree) {
       Tree *newtree = new Tree();
       newtree->type = VALUE;
diff --git a/unittest/commands/test_variables.cpp b/unittest/commands/test_variables.cpp
index eb533aee86..97f874a856 100644
--- a/unittest/commands/test_variables.cpp
+++ b/unittest/commands/test_variables.cpp
@@ -142,12 +142,13 @@ TEST_F(VariableTest, CreateDelete)
     command("variable ten2   uloop     4");
     command("variable ten3   uloop     4 pad");
     command("variable dummy  index     0");
+    command("variable file   equal     is_file(MYFILE)");
     END_HIDE_OUTPUT();
-    ASSERT_EQ(variable->nvar, 17);
+    ASSERT_EQ(variable->nvar, 18);
     BEGIN_HIDE_OUTPUT();
     command("variable dummy  delete");
     END_HIDE_OUTPUT();
-    ASSERT_EQ(variable->nvar, 16);
+    ASSERT_EQ(variable->nvar, 17);
     ASSERT_THAT(variable->retrieve("three"), StrEq("three"));
     variable->set_string("three", "four");
     ASSERT_THAT(variable->retrieve("three"), StrEq("four"));
@@ -158,6 +159,13 @@ TEST_F(VariableTest, CreateDelete)
     ASSERT_THAT(variable->retrieve("eight"), StrEq(""));
     variable->internal_set(variable->find("ten"), 2.5);
     ASSERT_THAT(variable->retrieve("ten"), StrEq("2.5"));
+    ASSERT_THAT(variable->retrieve("file"), StrEq("0"));
+    FILE *fp = fopen("MYFILE","w");
+    fputs(" ",fp);
+    fclose(fp);
+    ASSERT_THAT(variable->retrieve("file"), StrEq("1"));
+    unlink("MYFILE");
+    ASSERT_THAT(variable->retrieve("file"), StrEq("0"));
 
     ASSERT_EQ(variable->equalstyle(variable->find("one")), 0);
     ASSERT_EQ(variable->equalstyle(variable->find("two")), 1);
@@ -200,7 +208,7 @@ TEST_F(VariableTest, CreateDelete)
 
 TEST_F(VariableTest, AtomicSystem)
 {
-    command("atom_modify map array");
+    HIDE_OUTPUT([&] { command("atom_modify map array"); });
     atomic_system();
     file_vars();