Estimating VSHALE from GR remains the most popular standard technique to compute an initial curve of rock shaliness. However, the approach is far from perfect due to several criticisms:
 How reliable are the VSH estimates from GR if the rock matrix itself contains some even slightly radioactive minerals like feldspars?. What about if there are Kaolinite clay minerals with a low Gamma signal response?. If you are luckily enough, you may have a modern spectral Gamma Ray suite to apply a careful computation, better than using traditional total Gamma Ray.
 If the first guess linear Gamma Ray Index produces an overestimation of the VSH (very common), you may reduce it using well known corrections VSH=f(IGR), like Larionov for tertiary rocks or young rocks, Larionov for older rocks, or even Steiber and Clavier. But what is the best correction for your actual reservoir?
 What would happen if an increase in clay content is related to an increase in rock bulk density?
 To apply a Vshale estimate from GR, you must be sure you really have a clean sand and a real shale rock in the formation under study. But what would happen if just the whole clastic formation is shaly?, then you can not simply pick a value of GR_clean to input into the linear GR index equation. Likewise, are you sure do you have a real shale rock body in your formation?
If there are available modern logs, it is good to compare several Volume of Shale estimates. The Neutron Porosity to Density Porosity difference is perhaps after the GR, the next preferred method to estimate VSHALE.
Provided that the reservoir under study does not contain gas, shaly sands will show higher neutron porosity values than density porosity values. Conversely, if neutron porosity is smaller than density porosity (called gas crossover), it is a very strong indication of gas presence.
With no gas present in the formation, the neutron porosity to density porosity difference has a linear response with rock clay content. Even better, it does not depend upon GR, so it provides an independent reference estimator:
The GLS script program shown below to estimate VSH runs "out of the box". On a typical Windows installation just follow: [File] > [Open_GLS_Script] > [Open] > [UserHomeDir] > [GeolOil_DB] > [GLS_examples] > [neutr_densit_phi.gls] > [Run]. The script will create a new .LAS file with curves to compare VSH from GR, against VSH from neutrondensity porosity difference.
The figure below shows a GLS script program to estimate VSH using the neutrondensity method
The log curves plotted below show the result of the new generated .LAS file computed by the GLS script:

Track 1
Shows basic general and lithology information: The Gamma Ray curve GRGC (a rule of thumb of 65 API units threshold value was used to contrast in yellow and green, sandstones and shaly sandstones trends). The photoelectric factor PDPE, and the caliper CLDC. 
Track 2
Plots resistivity. A 17 Ohm.m threshold cutoff pink filling was used just as a preliminary guide to explore possible pay sand zones. 
Track 3
Highlights in light green the difference between neutron porosity and density porosity. The area highlighted itself is an excellent visual indicator of clay content. The thicker, the more volume of shale. The thinner, the cleaner rock. Furthermore, if there would be gas presence, the neutron to density crossing would happen clearly negative in the reverse direction. A bonus. 
Track 4
Compares three classic estimates for Vshale: The solid and thick green curve represents the standard linear Gamma Ray Index, it seems to have a trend to overestimate VSH in many cases. The dotted light green curve shows the Larionov youngshallow rocks VSH correction from IndexGR curve. It seems to subestimate the volume of shale. The solid red curve plots the VSH estimate using the neutron porosity to density porosity difference approach. It has a balanced and believable behaviour of rock shaliness, showing higher resolution details than the GR log.
The neutron to density difference method usually does not need corrections for matrix radioactivity, or VSH overestimation adjustments. It is indeed an excellent VSH indicator that should be computed any time the logs are available and the caliper shows a stable borehole, as the density porosity would be severely affected.
The figure below compares VSH computed as Linear GR index, Larionov and NeutronDensity difference.