diff --git a/include/particle_plt.hh b/include/particle_plt.hh
index 91cd4d9..81c5fcb 100644
--- a/include/particle_plt.hh
+++ b/include/particle_plt.hh
@@ -169,11 +169,16 @@ private:
                         D_xx_.kelem(i,j,k) *= norm;
                         D_yy_.kelem(i,j,k) *= norm;
                         D_zz_.kelem(i,j,k) *= norm;
+
+                        // spatially dependent correction to vfact = \dot{D_+}/D_+
+                        D_xy_.kelem(i,j,k) = 1.0/(0.25*(std::sqrt(1.+24*eval[2])-1.));
                     }
                 }
             }
         }
 
+        D_xy_.kelem(0,0,0) = 1.0;
+
     }
 
 public:
@@ -207,6 +212,12 @@ public:
         return D_zz_.get_cic_kspace({ix,iy,iz});
     }
 
+    inline ccomplex_t vfac_corr( std::array<size_t,3> ijk ) const
+    {
+        real_t ix = ijk[0]*mapratio_, iy = ijk[1]*mapratio_, iz = ijk[2]*mapratio_;
+        return D_xy_.get_cic_kspace({ix,iy,iz});
+    }
+
 };
 
 #if 0
diff --git a/src/ic_generator.cc b/src/ic_generator.cc
index 04682ad..90fbb0a 100644
--- a/src/ic_generator.cc
+++ b/src/ic_generator.cc
@@ -549,6 +549,9 @@ int Run( ConfigFile& the_config )
                                 tmp.kelem(idx) = vunit / boxlen * ( lg.gradient(idim,{i,j,k}) * phitot_v 
                                         + vfac3 * (lg.gradient(idimp,{i,j,k}) * A3[idimpp]->kelem(idx) - lg.gradient(idimpp,{i,j,k}) * A3[idimp]->kelem(idx)) );
 
+                                // correct velocity with PLT mode growth rate
+                                tmp.kelem(idx) *= lg.vfac_corr({i,j,k});
+
                                 if( bAddExternalTides ){
                                     // modify velocities with anisotropic expansion factor**2
                                     tmp.kelem(idx) *= std::pow(lss_aniso_alpha[idim],2.0);