diff --git a/README.md b/README.md
index 74e53c3a5c00b2cda703aa59e30b50c6b92956f3..8382022b0b0069248aa3dcae3d4cfb9732ddc9a4 100644
--- a/README.md
+++ b/README.md
@@ -200,6 +200,44 @@ python $codes/Becke_grid.py ${x}.scf.h5 chi.dat 4
 # Mainly for prepare density, gradient, laplacian. See [https://doi.org/10.1021/acs.jctc.0c00898] if you need to use it
 # We use it for the completeness.
 python $codes/rho_Becke_grid.py ${x}.scf.h5 Becke_grid.npy Becke_wgt.npy all 0
+
+
+examples/ch4_200_031/prep/mlcorr_prep.sh:
+
+#! /bin/bash
+
+# molecule name, here ch4
+x=$1
+
+# the folder of codes
+codes=../../../codes/eneden
+
+# the folder of prepared data
+data=../prep
+
+# generate Ajb（rg） matrix, see https://doi.org/10.1021/acs.jctc.0c00898
+# very time consuming, can be parallelize, or calculated on GPU (you need pycuda/skcuda)
+# 0,70000 means the first 70000 quadruture points for atom C,
+# 70000,140000 means the second 70000 quadruture points for atom H,
+# other H atoms are symmetric and are not needed to calculated
+# generate Ajb_200_031.C.npy Ajb_200_031.H.npy default use 8 cpu cores
+python $codes/prepare_Ajb.cpu.py $data/$x.scf.h5 $data/FCIDUMP_int2el.npy $data/$x.npz 0,70000 200_031.C.npy
+python $codes/prepare_Ajb.cpu.py $data/$x.scf.h5 $data/FCIDUMP_int2el.npy $data/$x.npz 70000,140000 200_031.H.npy
+
+# generate t(n) amplitudes for PT2, PT3, PT4
+# very time consuming, can be parallelize, or calculated on GPU (you need pytorch-gpu)
+# generate mp2_amp.npy mp3_amp.npy mp4_amp.npy (t(n)) default use 8 cpu cores
+python $codes/int2el_mp2-pytorch.cpu.py $data/$x.scf.h5 $data/FCIDUMP_int2el.npy
+python $codes/int2el_mp3-pytorch.cpu.py $data/$x.scf.h5 $data/FCIDUMP_int2el.npy
+python $codes/int2el_mp4-pytorch.cpu-j.mat.py $data/$x.scf.h5 $data/FCIDUMP_int2el.npy 
+
+# generate energy density for C and H
+for n in 2 3 4
+do
+  echo "Calculated PT$n (summing up of energy density for checking):"
+  python $codes/Ajb_mpn.pytorch.py $data/$x.scf.h5 mp${n}_amp.npy Ajb_200_031.C.npy $data/$x.npz 0,70000 200_031.C.mp${n}.npy
+  python $codes/Ajb_mpn.pytorch.py $data/$x.scf.h5 mp${n}_amp.npy Ajb_200_031.H.npy $data/$x.npz 70000,140000 200_031.H.mp${n}.npy
+done
 ```