anchor.diff:          n_53  = n_43*n_13
           n_83  = n_53*n
   
  -        x      = agr2/n_83
  -        dxdn   = -etthrd*x*inv_n
  -        dxdagr = 2.0d0*agr/n_83
  -        z      = tau/n_53 - cf
  -        dzdn   = -fvthrd*tau/n_83
  -        dzdtau = 1.0d0/n_53
  +        x       =  agr2/n_83
  +        dx_dn   = -etthrd*x*inv_n
  +        dx_dagr =  2.0d0*agr/n_83
  +
  +        z       =  tau/n_53 - cf
  +        dz_dn   = -fvthrd*tau/n_83
  +        dz_dtau =  1.0d0/n_53
           
  -********* VS98 Exchange *********
  +*       ***** VS98 Exchange *****
           aa    = -9.800683d-1
           bb    = -3.556788d-3
           cc    =  6.250326d-3
  
positive.code:         call nwpw_GVT4(aa,bb,cc,dd,ee,ff,alpha,alpha,
       >                 xrud,zrud,gammaud,GGud,dGud_dxud,dGud_dzud)
  
          eudc        = GGud*eud0c
          dGud_dnup   = dGud_dxud*dxup_dnup + dGud_dzud*dzup_dnup
          dGud_dagrup = dGud_dxud*dxup_dagrup 
          dGud_dtauup = dGud_dzud*dzup_dtauup
          dGud_dndn   = dGud_dxud*dxdn_dndn + dGud_dzud*dzdn_dndn
          dGud_dagrdn = dGud_dxud*dxdn_dagrdn 
          dGud_dtaudn = dGud_dzud*dzdn_dtaudn
  
          !fudc = n*GGud*eud0c
          fud_nupc    = GGud*eud0c 
       &                + nup*(dGud_dnup*eud0c + GGud*fud0_nupc)
          fud_dnupc   = nup*eud0c*dGud_dagrup
          fud_dtauupc = nup*eud0c*dGud_dtauup
          fud_ndnc    = GGud*eud0c 
       &                + ndn*(dGud_dndn*eud0c + GGud*fud0_ndnc)
          fud_dndnc   = ndn*eud0c*dGud_dagrdn
          fud_dtaudnc = ndn*eud0c*dGud_dtaudn
  
          ec       = euuc + eddc + eudc
          fnupc    = fuu_nupc + fud_nupc
          fdnupc   = fuu_dnupc + fud_dnupc
          fdtauupc = fuu_dtauupc + fud_dtauupc 
          fndnc    = fdd_ndnc + fud_ndnc
          fdndnc   = fdd_dndnc + fud_dndnc
          fdtaudnc = fdd_dtaudnc + fud_dtaudnc 
          
          xce(i)     = x_parameter*ex       + c_parameter*ec
          fn(i,1)    = x_parameter*fnupx    + c_parameter*fnupc
          fn(i,2)    = x_parameter*fndnx    + c_parameter*fndnc
  
negative.code:         ee    = -5.451559d-3
          ff    =  1.577575d-2
          alpha =  0.00515088d0
  *       ***** alpha-alpha *****
          gammaup = 1.0d0 + alpha*(xup + zup)
          xrup    = xup/gammaup
          zrup    = zup/gammaup
  
          call nwpw_vs98c_ss(aa,bb,cc,dd,ee,ff,alpha,nup,agrup,tauup,
       >               xup,zup,gammaup,xrup,zrup,
       >               dxup_dnup,dxup_dagrup,dzup_dnup,dzup_dtauup,
       >               euu0c,deuu0c_dnup,
       >               euuc,dfuuc_dnup,dfuuc_dagrup,dfuuc_dtauup)
  
  *       ***** beta-beta *****
          gammadn = 1.0d0 + alpha*(xdn + zdn)
          xrdn    = xdn/gammadn
          zrdn    = zdn/gammadn
  
          call nwpw_vs98c_ss(aa,bb,cc,dd,ee,ff,alpha,ndn,agrdn,taudn,
       >               xdn,zdn,gammadn,xrdn,zrdn,
       >               dxdn_dndn,dxdn_dagrdn,dzdn_dndn,dzdn_dtaudn,
       >               edd0c,dedd0c_dndn,
       >               eddc,dfddc_dndn,dfddc_dagrdn,dfddc_dtaudn)
  
  *       ***** opposite-spin *****
          aa    =  7.035010d-1
          bb    =  7.694574d-3
          cc    =  5.152765d-2
  